################################################################################
#
# Triaxial test. Axial strain rate is prescribed and transverse prestress.
# Test is possible on prism or cylinder
# An independent c++ engine may be created from this script in the future.
#
################################################################################
import os
import sys
#import cpler2 (mpi binding)
from libyade import *
#cpler2.cpl2.init() (init comm mpi)
class Cplsim():
  def __init__(self):

    concMat = O.materials.append(CpmMat(
        young=20e9,frictionAngle=1.2,poisson=0.24,sigmaT=1.5e6,
        epsCrackOnset=1e-04,relDuctility=30
    ))
    frictMat = O.materials.append(FrictMat(
        young=20e9,poisson=0.24,frictionAngle=1.2
    ))

    # spheres
    height = 5e-3; width = 2e-3; rParticle = 0.075e-3; nw = 24; nh=15; bcCoeff=5; strainRate = -100
    pred = pack.inCylinder((0,0,0),(0,0,height),.5*width)
    sp=SpherePack()
    sp = pack.randomDensePack(pred,spheresInCell=2000,radius=rParticle,memoizeDb='/tmp/triaxTestOnCylinder.sqlite',returnSpherePack=True)
    spheres=sp.toSimulation(color=(0,1,1),material=concMat)
    self.preStress = -3e6
    # bottom and top of specimen. Will have prescribed velocity
    bot = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]<rParticle*bcCoeff]
    top = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
    vel = strainRate*(height-rParticle*2*bcCoeff)
    for s in bot:
        s.shape.color = (1,0,0)
        s.state.blockedDOFs = 'xyzXYZ'
        s.state.vel = (0,0,-vel)
    for s in top:
        s.shape.color = Vector3(0,1,0)
        s.state.blockedDOFs = 'xyzXYZ'
        s.state.vel = (0,0,vel)

    # facets
    self.facets = []
    rCyl2 = .5*width / cos(pi/float(nw))
    for r in xrange(nw):
        for h in xrange(nh):
                v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
                v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
                v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
                v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
                f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)
                f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)
                self.facets.extend((f1,f2))

    O.bodies.append(self.facets)
    mass = O.bodies[0].state.mass
    for f in self.facets:
        f.state.mass = mass
        f.state.blockedDOFs = 'XYZz'


    O.dt=.5*utils.PWaveTimeStep()
    enlargeFactor=1.5
    O.engines=[
        ForceResetter(),
        InsertionSortCollider([
            Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
            Bo1_Facet_Aabb()
        ]),
        InteractionLoop(
            [
                Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ss2d3dg'),
                Ig2_Facet_Sphere_ScGeom(),
            ],
            [
                Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=O.iter+5),
                Ip2_FrictMat_CpmMat_FrictPhys(),
                Ip2_FrictMat_FrictMat_FrictPhys(),
            ],
            [
                Law2_ScGeom_CpmPhys_Cpm(),
                Law2_ScGeom_FrictPhys_CundallStrack(),
            ],
        ),
        PyRunner(iterPeriod=1,command="sim1.addForces()"),
        NewtonIntegrator(damping=.3),
        CpmStateUpdater(iterPeriod=50,label='cpmStateUpdater'),
        #PyRunner(command='plotAddData()',iterPeriod=10),
        #PyRunner(iterPeriod=50,command='stopIfDamaged()'),
    ]

  # run one step
    O.step()

  # reset interaction detection enlargement
    bo1s.aabbEnlargeFactor=ss2d3dg.interactionDetectionFactor=1.0


  def addForces(self):
        print "stp = ", O.iter
    for f in O.bodies:
        #n = f.shape.normal
        #a = f.shape.area
                if type(f.shape)==Sphere:
                 fluvel = Vector3(0.001,0.0,0.0) #some dummy variables to show performance loss
                 bvel = f.state.vel#some dummy variables to show performance loss
                 flurot = Vector3(0,0,0.005)#some dummy variables to show performance loss
                 fx = 6*3.14*f.shape.radius*0.001*(fluvel-bvel) #some dummy variables to show performance loss
                 mx = 8*3.14*(f.shape.radius**3)*(flurot-f.state.angVel)  #some dummy variables to show performance loss
                 ### addforce
         O.forces.addF(f.id,fx)  #if you add a constant force things run smmothly with max performance.
                 O.forces.addT(f.id,mx)



  def irun(self,num):
    print O.time
    O.run(num,1)


if __name__ == "__main__":
  sim1 = Cplsim()
  sim1.irun(5000)