model new
model title 'Triaxial Test' 
model random 10001
model domain extent -2 2 -2 2 -4 4 condition destroy
model large-strain on
model mechanical timestep scale

[rad = 1]
[height = 4]
[segments = 6]
[halfLen = height/2.0]
[freeRegion = (height/2.0)*0.8]

geometry edge create by-positions (-1 -1 -2) (-1 1 -2) (1 1 -2) (1 -1 -2) (-1 -1 -2)

geometry generate from-edges extrude (0,0,[height]) segments [segments*2]

structure shell import from-geometry 'Default' element-type dkt-cst 
structure node group 'middle' range position-z [-freeRegion] [freeRegion]
structure node group 'top' range position-z [freeRegion] [freeRegion+0.4]
structure node group 'bot' range position-z [-freeRegion-0.4] [-freeRegion]
structure shell group 'middle' range position-z [-freeRegion] [freeRegion]
structure shell property isotropic (2e7, 0.2) thick 0.04 

model cycle 0
;In order to use the STRUCTURE APPLY command 
;set the local system of each structural
;element to be pointing to the center of the triaxial cell
fish define setLocalSystem
    loop foreach local s struct.node.list()
        local p = struct.node.pos(s)
        local nid = struct.node.id.component(s)
        local mvec = vector(0,0,comp.z(p))
        zdir = math.unit(p-mvec)
        ydir = vector(0,0,1)
        command
            structure node system-local z [zdir] y [ydir] ...
                      range component-id [nid]
        endcommand
    endloop
   command
     structure node fix system-local
   endcommand
end
[setLocalSystem]
structure damp local
structure node fix velocity rotation

wall-structure create

[rad2 = rad*1.3]

wall generate name 'platenTop' polygon ([-rad2],[-rad2],[halfLen]) ...
                                       ([rad2],[-rad2],[halfLen])  ...
                                       ([rad2],[rad2],[halfLen])   ...
                                       ([-rad2],[rad2],[halfLen])
wall generate name 'platenBottom' polygon ([-rad2],[-rad2],[-halfLen]) ...
                                          ([rad2],[-rad2],[-halfLen])  ...
                                          ([rad2],[rad2],[-halfLen])   ...
                                          ([-rad2],[rad2],[-halfLen])                                          
wall resolution full
wall attribute cutoff-angle 20
model save 'first'
model restore 'first'

ball distribute box [-rad] [rad] [-rad] [rad] [-halfLen] [halfLen] ...
                porosity 0.3684 bin 1 radius 0.043 0.175
ball attribute density 2600 damp 0.3

contact cmat default model linear ...
                 property kn 9.6e8 ks 3.84e7
contact cmat default type ball-facet model linear ...
                  property kn 1e7 ks 1e7
model cycle 1000 calm 100

contact cmat default model rrlinear ...
                   property kn 9.6e8 ks 3.84e7 fric 0.36 ...
                   rr_fric 0.1
contact cmat default type ball-facet model linear ...
                  property kn 1e7 ks 1e7

contact cmat apply
model solve ratio-average 1e-5

model calm
ball fix velocity spin
model cycle 2
ball free velocity spin

structure node free velocity rotation range group 'middle'

[platenTop = wall.find('platenTop')]
[platenBottom = wall.find('platenBottom')]
[failureStress = 0]
[currentStress = 0]
[failureStrain = 0]
[area = rad*rad*2*2]
fish define stress
    local topForce = math.abs(comp.z(wall.force.contact(platenTop))) 
    local botForce = math.abs(comp.z(wall.force.contact(platenBottom)))
    currentStress = 0.5*(topForce+botForce)/area
    stress = currentStress
    strain = (height - (comp.z(wall.pos(platenTop)) - ...
              comp.z(wall.pos(platenBottom))))/height * 100
    if failureStress <= currentStress
        failureStress = currentStress
        failureStrain = strain
    endif
end

fish define halt
    halt = 0
    if currentStress < failureStress * 0.85
        halt = 1
    endif
end

fish define rampUp(beginIn,ending,increment)
    command
        ball attribute displacement (0,0,0)
        structure node initialize displacement (0,0,0)
    endcommand
    begin = beginIn
    loop while (math.abs(begin) < math.abs(ending))
        begin = begin + increment
        command
            structure shell apply [begin] range group 'middle'
            wall servo force (0,0,[begin*area]) activate true ...
                 range name 'platenTop'
            wall servo force (0,0,[-begin*area]) activate true ...
                 range name 'platenBottom'
model cycle 200
            model calm
        endcommand
    endloop
    command

model cycle 1000
        wall servo activate false
        wall attribute velocity (0,0,0) range name 'platenTop'
        wall attribute velocity (0,0,0) range name 'platenBottom'
    endcommand
end
model save '11'
model restore '11'
[platenVel = 0.000003]

model save 'beforeApplication'
model restore 'beforeApplication'
[rampUp(0,-1e5,-1e3)]
model save 'to1e4'
model restore 'to1e4'
wall attribute velocity-z [-platenVel] range name 'platenTop'
wall attribute velocity-z [platenVel] range name 'platenBottom'
ball attribute displacement (0,0,0)
structure node initialize displacement (0,0,0)
fish history stress
fish history strain
model solve fish-halt halt
model save 'triaxial1e4'
[io.out(string(failureStress) + 'Pa ')]
[io.out('at' + string(failureStrain) + '% strain')]