plastic energy dissipation

### this simulation models impact between a boulder and a 2D soil layer.
### the soil layer is divided into several parts and plastic energy dissipation inside each parts are tracked in order to investigate the impact mechanisms.



##### part 1 : generate a soil sample by gravity deposit

from yade import pack,plot
from yade import utils
import random


finalFricDegree = 30
O.trackEnergy=True
young=100e6
rmean=0.1

rmin=0.05
rmax=0.15

mn=Vector3(0,0,-2*rmean)
mx=Vector3(30*2*rmean,30*2*rmean,2*rmean)

number=500
fallingHeight=5

O.materials.append(CohFrictMat(young=100e6,poisson=0.3,frictionAngle=radians(0),density=2600,normalCohesion=0,shearCohesion=0,alphaKr=0.5,alphaKtw=0.01,momentRotationLaw=False,isCohesive=False,etaRoll=-1,label='spheres'))

O.materials.append(FrictMat(young=100e6,poisson=0.3,frictionAngle=0,density=0,label='walls'))


walls=utils.aabbWalls([mn,mx],thickness=0.5*rmean, oversizeFactor=1.0, material='walls')
wallIds=O.bodies.append(walls)



ssp=pack.SpherePack()

for i in range(1,number):
    r =rmin + (rmax-rmin)*(-0.15*log(1-random.random()))**1.33333
    if r>rmax:
     r =rmin + (rmax-rmin)*(-0.15*log(1-random.random()))**1.33333
    mbegin=(0,0,-1*r)
    mend=(mx[0],mx[1],1*r)
    ssp.makeCloud(minCorner=mbegin,maxCorner=mend,rMean=r,rRelFuzz=0,num=1)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in ssp])


newton=NewtonIntegrator(damping=0.9,kinSplit=True,gravity=(0,0,0))

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack(label='law1'),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
                              useIncrementalForm=False,
                              always_use_moment_law=False,
                              label='law2')]
    ),
    GravityEngine(gravity=(0,-9.81,0)),
    newton

]

for b in O.bodies:
    if isinstance(b.shape,Sphere): b.state.blockedDOFs='zXY' # blocked movement along Z


O.dt=utils.PWaveTimeStep()

O.run(150000, True)

## reset boundaries
Mz=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
h=Mz/4

mn=(0,0,-2*rmean)
mxafter=Vector3(mx[0],Mz+2*rmean,2*rmean)

O.bodies.erase(0)
O.bodies.erase(1)
O.bodies.erase(2)
O.bodies.erase(3)
O.bodies.erase(4)
O.bodies.erase(5)

O.bodies.append(utils.aabbWalls([mn,mxafter],thickness=0.1*rmean, oversizeFactor=1.0, material='walls'))

O.bodies.erase(3)
O.bodies.erase(4)
O.bodies.erase(5)

O.run(10000, True)

## add friction and moment law
O.materials[0].frictionAngle=radians(finalFricDegree)
O.materials[0].momentRotationLaw=True
O.engines[2].lawDispatcher.functors[1].always_use_moment_law = True

O.engines[4].damping=0.000

O.run(80000, True)
O.saveTmp()

##### part 2: divide the soil layer into different parts


mxafter=Vector3(mx[0],Mz,2*rmean)

centerx=0.5*mxafter[0]
centery=mxafter[1]



O.materials.append(CohFrictMat(young=100e6,poisson=0.3,frictionAngle=radians(finalFricDegree),density=2600,normalCohesion=0,shearCohesion=0,alphaKr=0.5,alphaKtw=0.01,
        momentRotationLaw=True,isCohesive=False,etaRoll=-1,label='heart'))

for i in range(1,4):
    for j in range(1,3):
        O.materials.append(CohFrictMat(young=100e6,poisson=0.3,frictionAngle=radians(finalFricDegree),density=2600,normalCohesion=0,shearCohesion=0,alphaKr=0.5,alphaKtw=0.01,
        momentRotationLaw=True,isCohesive=False,etaRoll=-1,label='crown_'+str(i)+str(j)))


import numpy,random
a=numpy.zeros(6).reshape((3, 2))


for i in O.bodies:
  x=i.state.pos[0]
  y=i.state.pos[1]
  D=sqrt((x-centerx)**2+(y-centery)**2)
  xy=atan2((y-centery),(x-centerx))
  k=0
  if (D/h <= 1):
      i.material=O.materials['heart']
      i.shape.color=Vector3(1,0,0)
  if (1 < D/h < 4):
       j=int(D/h)
       if -2*pi/3 <= xy <= -pi/3 :
           k=1
       else:
           k=0
       a[j-1][k]+=1
       i.material=O.materials['crown_'+str(j)+str(k+1)]


## mark the soil layers with different colors
color=Vector3(random.random(),random.random(),random.random())

for i in O.bodies:
    i.shape.color=color

for i in O.bodies:

    if i.material==O.materials['crown_11']:
       i.shape.color=Vector3(0.2,0.5,0)

    if i.material==O.materials['crown_12']:
       i.shape.color=Vector3(0.6,0.2,1)

    if i.material==O.materials['crown_21']:
       i.shape.color=Vector3(1,0.5,0.6)

    if i.material==O.materials['crown_22']:
       i.shape.color=Vector3(1,1,0)

    if i.material==O.materials['crown_31']:
       i.shape.color=Vector3(0.2,0.9,0.3)

    if i.material==O.materials['crown_32']:

       i.shape.color=Vector3(0.6,1,1)


    if i.material==O.materials['heart']:
       i.shape.color=Vector3(1,0,0)


Gl1_Sphere.stripes=1
O.run(10000,True)
O.save('2dlayer.xml')

#### part 3: add the boulder and conduct impact modelling, track energy dissipations inside each soil layer

boulderradius=3*rmean

sp2=pack.SpherePack()

mbegin=(0.5*mxafter[0]*1.0-boulderradius,mxafter[1],-1*boulderradius)
mend=(0.5*mxafter[0]*1.0+boulderradius,mxafter[1]+2
*boulderradius,1*boulderradius)
sp2.makeCloud(minCorner=mbegin,maxCorner=mend,rMean=boulderradius,rRelFuzz=0,num=1)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp2])

O.bodies[3].state.blockedDOFs='zXY'

a=sin(0*3.14/180)
b=cos(0*3.14/180)
O.bodies[3].state.vel=Vector3(-sqrt(2*9.8*fallingHeight)*a,-sqrt(2*9.8*fallingHeight)*b,0)

O.dt=0.1*utils.PWaveTimeStep()

for i in range (1,200):
    O.run(100,True)
    if -O.forces.f(3)[1]< O.bodies[3].state.mass*9.81:
      tt=O.time
      break

## time when the impact begins
tt=O.time

O.saveTmp()

aa=dict(O.energy.items())
E=aa['plastDissip']
f=open('eptrack.txt','w')
ff=open('EkEstrack.txt','w')
p_11=0
p_12=0
p_21=0
p_22=0
p_31=0
p_32=0
p=0
psum=0

for k in range (1,4000):

    O.run(1,True)
    nEs=0
    sEs=0
    Emgh=0
    Ek=0
    aa=dict(O.energy.items())
    v=O.bodies[3].state.vel.norm()

    for i in O.bodies:
            if isinstance(i.shape, Sphere):
                   Emgh=Emgh+(i.state.mass)*(i.state.pos[1])*9.81
                   Ek=Ek+0.5*i.state.mass*i.state.vel.norm()**2+0.5*i.state.inertia[0]*i.state.angVel.norm()**2

    for i in O.interactions:

    Ft=i.phys.shearForce
    du=i.geom.shearInc
    kt=i.phys.ks
    kn=i.phys.kn
    Fn=i.phys.normalForce
    nEs=nEs+0.5*Fn.squaredNorm()/kn
    sEs=sEs+0.5*Ft.squaredNorm()/kt

    Ftt=Ft-kt*du
    maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
    if(Ftt.norm() > maxFs):
            ratio = maxFs / Ftt.norm()
            trialForce=Ftt
            Ftt *= ratio
            e=((1/kt)*(trialForce-Ftt)).dot(Ftt)
            psum+=e
            if (O.bodies[i.id1].material== O.materials['crown_11']) or (O.bodies[i.id2].material== O.materials['crown_11']):
                    p_11+=e
                elif (O.bodies[i.id1].material== O.materials['crown_12']) or (O.bodies[i.id2].material== O.materials['crown_12']):
                    p_12+=e
                elif (O.bodies[i.id1].material== O.materials['crown_21']) or (O.bodies[i.id2].material== O.materials['crown_21']):
                    p_21+=e
                elif (O.bodies[i.id1].material== O.materials['crown_22']) or (O.bodies[i.id2].material== O.materials['crown_22']):
                    p_22+=e
                elif (O.bodies[i.id1].material== O.materials['crown_31']) or (O.bodies[i.id2].material== O.materials['crown_31']):
                    p_31+=e
                elif (O.bodies[i.id1].material== O.materials['crown_32']) or (O.bodies[i.id2].material== O.materials['crown_32']):
                    p_32+=e
                else:
                    p+=e
    f.write('%f %f %f %f %f %f %f %f %f %f\n' %(O.time-tt,p_11,p_12,p_21,p_22,p_31,p_32,p,psum,aa['plastDissip']-E))
    ff.write('%f %f %f %f %f %f %f %f %f \n' %(O.time-tt,v,Emgh,Ek,nEs,sEs,utils.kineticEnergy(),law2.normElastEnergy(),law2.shearElastEnergy()))

f.close()
ff.close()