##Import Libraries

from yade import pack,plot,qt
import math
import random as rand
import numpy as np


servermode=False
if servermode:
    qt.View()
O.dt=1e-6

##############################
## PARAMETERS OF SIMULATION ##
##############################



### Configuration: inclined channel
slope = math.radians(31)    #Slope of the channe, in radian
length = 0.3    #Length of the channe, in m
width = 0.01    #Width of the channel, in m
height = 0.1    #Height of the channel, in m
gravityVector = Vector3(9.81*sin(slope),0.0,-9.81*cos(slope))   #Gravity vector to consider a channel inclined with slope angle 'slope'

### Materials
# Channel

#Fluid
fluidDensity= 1200. #Density of the fluid, in kg/m3

# Debris
diameterPart = 0.001    #Diameter of the particles, in m
O.materials.append(
        FrictViscoMat(
            betan=0.5,  #(not defined yet)
            density=1800.,  #entre 1800 et 2300 si densite apparente
            frictionAngle=math.radians(35.),
            poisson=.5,#Ratio between kn and ks (not defined yet)
            young=1e8,#Not the real young but contact stiffnesse kn (not defined yet)
            label='DebrisMat'
        ))

########################
## FRAMEWORK CREATION ##
########################
#
# Channel creation
O.bodies.append(
    geom.facetBox(
    (length/2.,width/2.,height/2.),
    (length/2.,width/2.,height/2.),
    Quaternion((1,0,0),0),
    wallMask=63,
    material='DebrisMat'))
O.bodies.append(
    geom.facetBox(
    (0,width/2.,height/2.),
    (3./100,width/2.,height/2.),
    Quaternion((1,0,0),0),
    wallMask=2,
    material='DebrisMat'))

rigWall=O.bodies.append(
        aabbWalls(
        [(0.9*length,0,0),(0.9*length,width,0.015)],
        thickness=0.001,
        material='DebrisMat',
        oversizeFactor=1.0,
        color=(204, 0, 0)))
for i in rigWall[1:]:
        O.bodies.erase(i)

# Creation of a pack of sphere
setclump=False

partCloud = pack.SpherePack()
partNumber = 5000
partCloud.makeCloud(minCorner=(0.01,0,0),maxCorner=(3./100,width,0.01),rRelFuzz=diameterPart/10, rMean=diameterPart/2.0)
spIds=partCloud.toSimulation(material='DebrisMat',color=[205./255,175./255,149./255])

if setclump:
    templates= []
    relRadList1=[diameterPart/2,diameterPart/2]
    relPosList1=[[diameterPart/2,0,0],[diameterPart,0,0]]
    templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
    ClIds=O.bodies.replaceByClumps(templates,[1])


#####################
## SIMULATION LOOP ##
#####################
Factor=2.  #(test value)
O.engines=[
    ForceResetter(),
    InsertionSortCollider(
    [
    Bo1_Sphere_Aabb(aabbEnlargeFactor=Factor),
    Bo1_Facet_Aabb(),
    Bo1_Box_Aabb(),
    Bo1_Wall_Aabb()
    ],
    label="ISC"
    ),
    InteractionLoop(
    # Geometries
    [
    Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=Factor),
    Ig2_Facet_Sphere_ScGeom(),
    Ig2_Box_Sphere_ScGeom()
    ],

    # Physics
    [
#   Ip2_FrictMat_FrictMat_FrictPhys()
    Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys()
#   Ip2_FrictMat_FrictViscoMat_FrictViscoPhys()
    ],

    # Laws
    [
    Law2_ScGeom_FrictViscoPhys_CundallStrackVisco(),
    ]
    ),

#   GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    NewtonIntegrator(gravity=gravityVector, damping = 0., label='newton'),
    VTKRecorder(fileName='Results/test/test-',iterPeriod=10000),
    # Launching Pyhton functions inside the loop
    PyRunner(command='stopme()', iterPeriod=100),
#   PyRunner(command='print(unbalancedForce())',iterPeriod=500),
    PyRunner(command='addPlotData()',iterPeriod=500),
    PyRunner(command='gravityDeposit()',iterPeriod=500,label='gravDep',dead=False),
    ]


#########################
## FUNCTION DEFINITION ##
#########################

def gravityDeposit():
   if (O.iter>100000 and kineticEnergy()<.05):
    if servermode:
        O.pause()
    O.bodies.erase(12)
    O.bodies.erase(13)
    gravDep.dead=True
    newton.damping=0
    print "Gravity deposit finished & gate opened"

#and unbalancedForce()<.06 and O.energy.items()[0][1]<0.005
def checkUnbalanced():
   if (max([b.state.vel.norm() for b in O.bodies])<0.01 and unbalancedForce()<.05):
    O.pause();plot.saveDataTxt('bbb1.txt.bz2');
    print "stab cond satisf"

def addPlotData():
   plot.addData(i=O.iter, unbalanced=unbalancedForce(), **O.energy)


O.trackEnergy=True
#O.dt=PWaveTimeStep()

plot.plots={'i':('unbalanced',None,O.energy.keys)}
plot.plot()


O.saveTmp()