Ep code

E=O.energy['plastDissip']
e=0
energy=open('energy.txt','w')
for k in range (0,100):
    a=open('state_1.txt','w')
    b=open('state_2.txt','w')
    for i in O.interactions:
          Ft=i.phys.shearForce
          a.write('%d %d %f %f %f\n' %(i.id1,i.id2,Ft[0],Ft[1],Ft[2]))

    a.close()

    O.step()

    for i in O.interactions:
        Ftt=i.phys.shearForce
        du=i.geom.shearInc
        kt=i.phys.ks
        Fn=i.phys.normalForce
        maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
        b.write('%d %d %f %f %f %f %f %f %f %f\n' %(i.id1,i.id2,Ftt[0],Ftt[1],Ftt[2],maxFs,du[0],du[1],du[2],kt))


    b.close()

    x = numpy.loadtxt('state_1.txt')
    y = numpy.loadtxt('state_2.txt')

    for i in range (0,len(y)):
        Ftt=Vector3(y[i][2],y[i][3],y[i][4])
        maxFs=y[i][5]
        du=Vector3(y[i][6],y[i][7],y[i][8])
        kt=y[i][9]
        Ft=Vector3(0,0,0)
        for j in range (0,len(x)):
            if y[i][0]==x[j][0] and y[i][1]==x[j][1]:
                Ft=Vector3(x[j][2],x[j][3],x[j][4])
        Ftc=Ft-kt*du
        if Ftc.norm()> maxFs:
                #Ftt=Ftt/Ftt.norm()*maxFs
                due=-(Ftt-Ft)/kt
                #e+=Ftt.norm()*(du.norm()-due.norm())
                e+=Ftt.norm()*(du-due).norm()
    energy.write('%f %f %f\n'%(O.time,e,O.energy['plastDissip']-E))


energy.close()