Yade

# -*- coding: utf-8
from yade import qt,plot
from yade.gridpfacet import *
import math
math.pi

##################
### PARAMETERS ###
##################
phi=20.
E=5.*1e10
color=[255./255.,102./255.,0./255.]
r=0.05
density_sp = 2700

# position of imported mesh
xpafet=7.5
ypafet=1.1
zpafet=7.5

################
### ENGINES  ###
################
O.engines=[
    ForceResetter(),
    InsertionSortCollider([
        Bo1_PFacet_Aabb(),
        Bo1_Sphere_Aabb(),
    ]),
    InteractionLoop([
        Ig2_GridNode_GridNode_GridNodeGeom6D(),
        Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
#       Ig2_Sphere_GridConnection_ScGridCoGeom(),
        Ig2_GridConnection_PFacet_ScGeom(),
        Ig2_Sphere_PFacet_ScGridCoGeom(),
        Ig2_PFacet_PFacet_ScGeom(),
    ],
    [
        Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True),
        Ip2_FrictMat_FrictMat_FrictPhys()
    ],
    [
        Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),
        Law2_ScGeom_FrictPhys_CundallStrack(),
        Law2_ScGridCoGeom_FrictPhys_CundallStrack(),
        Law2_GridCoGridCoGeom_FrictPhys_CundallStrack()
    ]
    ),
    GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.5,label='ts'),
    NewtonIntegrator(gravity=(0,-9.81,0),damping=0.0,label='newton'),
    PyRunner(iterPeriod=200,command='history()'),
]

################
### MATERIAL ###
################
O.materials.append( CohFrictMat( young=E/100,poisson=0.3,density=2700,frictionAngle=radians(phi),normalCohesion=1e100,shearCohesion=1e100,momentRotationLaw=True,label='gridNodeMat1' ) )  # material to create the gridConnections
O.materials.append( CohFrictMat( young=E/100,poisson=0.3,density=1,frictionAngle=radians(phi),normalCohesion=1e100,shearCohesion=1e100,momentRotationLaw=True,label='gridConMat1' ) )
O.materials.append( CohFrictMat( young=E,poisson=0.3,density=2700,frictionAngle=radians(phi),normalCohesion=1e100,shearCohesion=1e100,momentRotationLaw=True,label='gridNodeMat2' ) )
O.materials.append( CohFrictMat( young=E,poisson=0.3,density=2700,frictionAngle=radians(phi),normalCohesion=1e100,shearCohesion=1e100,momentRotationLaw=True,label='gridConMat2' ) )
O.materials.append( FrictMat( young=E,poisson=0.3,density=2700,frictionAngle=radians(phi),label='pFacetMat') )  # material for general interactions
O.materials.append( FrictMat( young=E,poisson=0.3,density=density_sp,frictionAngle=radians(phi),label='sphereMat' ) )

########################################
### GENERATE THE WALL AND THE SPHERE ###
########################################

# FIXED WALL
N=30
IdNodes = []
for i in range(N):
    for j in range(N):
        IdNodes.append(O.bodies.append( gridNode([5*i*1e-1,0,5*j*1e-1],r,wire=False,fixed=False,material='gridNodeMat1',color=color) ))
        if i == 0 or i == (N-1) or j == 0 or j == (N-1):
            O.bodies[IdNodes[-1]].state.blockedDOFs='xyzXYZ'
for i in range(N*N-1):
    if i<(N*N-N):
        IdNodes.append(O.bodies.append( gridConnection(IdNodes[i],IdNodes[i+N],r,color=color,material='gridConMat1')))
    if (i+1)%N==0:
        continue
    IdNodes.append(O.bodies.append( gridConnection(IdNodes[i],IdNodes[i+1],r,color=color,material='gridConMat1') ))
    if i<(N*N-N):
        IdNodes.append(O.bodies.append( gridConnection(IdNodes[i],IdNodes[i+1+N],r,color=color,material='gridConMat1')))
for i in range(N*N-1):
    if (i+1)%N==0:
        continue
    if i<(N*N-N-1):
        IdNodes.append(O.bodies.append( pfacet(IdNodes[i],IdNodes[i+1],IdNodes[i+N+1],wire=False,material='pFacetMat',color=color) ))
        IdNodes.append(O.bodies.append( pfacet(IdNodes[i],IdNodes[i+N],IdNodes[i+1+N],wire=False,material='pFacetMat',color=color) ))

# LOWER MESH
N=30
IdNodes2 = []
for i in range(N):
    for j in range(N):
        IdNodes2.append(O.bodies.append( gridNode([5*i*1e-1,-1,5*j*1e-1],r,wire=False,fixed=False,material='gridNodeMat2',color=color) ))
        O.bodies[IdNodes2[-1]].state.blockedDOFs='xyzXYZ'
for i in range(N*N):
    m = O.bodies[IdNodes[i]].id
    n = O.bodies[IdNodes2[i]].id
    print (m,n)
    O.bodies.append( gridConnection(IdNodes[i],IdNodes2[i],r,color=color,material='gridConMat2') )


O.bodies.append(sphere(Vector3(xpafet,ypafet,zpafet), radius=20*r, wire=False, fixed=False, material='sphereMat', color=[1,0,0]))

############
### PLOT ###
############
init_pos = 0 #change it
m = O.bodies[-1].state.mass
#m = density_sp * 4/3 * pi * (20*r)**(3)
def history():
  xyz=[]
  pot = 0
  pot+=m * (9.81) * (O.bodies[-1].state.pos[1]-ypafet)
  kin = 0
  vel = O.bodies[-1].state.vel.norm()
  kin+=(m/2) * vel**(2)
  Force = 0

  for i in range(5163): #impatto in prossimità di un nodo
    d = dist(O.bodies[IdNodes[i]],O.bodies[-1])
    if d <= (21.3)*r: #check if the value can be lower
      print (i)
      break

  for k in [0,1,2]:
    ksum=0
    ksum+=O.bodies[-1].state.pos[k]
    xyz.append(ksum)  # take average value as reference
  plot.addData(i=O.iter,t=O.time,time=O.time,tt=O.time,x=xyz[0],y=xyz[1],z=xyz[2],pot_en=pot,kin_en=kin,tot_en=pot+kin,pen=abs(O.bodies[IdNodes[i]].state.pos[1]-init_pos),Force=O.forces.f(O.bodies[IdNodes[i]].id)[1])
plot.plots={'x':'y','t':('pot_en','kin_en','tot_en'),'time':'pen','tt':'Force'}
plot.plot(subPlots = False)

def dist(a,b):
    return math.sqrt((a.state.pos[0]-b.state.pos[0])**2+(a.state.pos[1]-b.state.pos[1])**2+(a.state.pos[2]-b.state.pos[2])**2)


##########
## VIEW ##
##########

qt.Controller()
qtv = qt.View()
qtr = qt.Renderer()
qtr.light2=True
qtr.lightPos=Vector3(1200,1500,500)
qtr.bgColor=[1,1,1]

O.saveTmp()