;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;                                       UCSD 2D Shear Beam Example                                       ;
;                                             March 08, 2017                                             ;
;--------------------------------------------------------------------------------------------------------;
set processors 1
;------------MODEL OPTIONS-------------------------------------------------------------------------------;
config dynamic fluid cppudm
gen zone brick &
p0 0 0 0 p1 1 0 0 p2 0 1 0 p3 0 0 10 &
size 1,1,10

model fluid fl_iso rang z = 0,10

set gravity 0 0 -9.81
set fluid off
set dyn off

ini fmod 2.2e6
ini fdensity 1.0

;------------MATERIAL PROPERTIES-------------------------------------------------------------------------;
;------------ASSIGN UCSD DEMO MODEL----------------------------------------------------------------------;
model ucsdsand3demo range name=Layer_4 z=0,10
prop density 1.99 refShearModul 7.3700E4 refBulkModul 1.96800E5 frictionAng 3.03E1  peakShearStra 0.1 refPress 1.01E2 pressDependCoe 5.E-1 PTAng 2.53E1 contrac1 0.035 contrac2 3 contrac3 0.2 dilat1 0.15 dilat2 3 dilat3 0.2 noYieldSurf 20 grTab 0 liquefac1 1 liquefac2 0 pAtm 1.01E2 cohesi 1.73E0 mType 0 b1 5.E-1 b2 3.0E1 b3 0.E0 b4 1.E-3 range=nrange Layer_4

prop Perm 6.6E-5 poro 5.E-1 range=nrange Layer_4

initial sxx -178.184385279 rang z=0,1
initial sxx -159.428134197 rang z=1,2
initial sxx -140.671883115 rang z=2,3
initial sxx -121.915632033 rang z=3,4
initial sxx -103.159380951 rang z=4,5
initial sxx -84.403129869 rang z=5,6
initial sxx -65.646878787 rang z=6,7
initial sxx -46.890627705 rang z=7,8
initial sxx -28.134376623 rang z=8,9
initial sxx -9.378125541 rang z=9,10

initial syy -178.184385279 rang z=0,1
initial syy -159.428134197 rang z=1,2
initial syy -140.671883115 rang z=2,3
initial syy -121.915632033 rang z=3,4
initial syy -103.159380951 rang z=4,5
initial syy -84.403129869 rang z=5,6
initial syy -65.646878787 rang z=6,7
initial syy -46.890627705 rang z=7,8
initial syy -28.134376623 rang z=8,9
initial syy -9.378125541 rang z=9,10

initial szz -185.45805 rang z=0,1
initial szz -165.93615 rang z=1,2
initial szz -146.41425 rang z=2,3
initial szz -126.89235 rang z=3,4
initial szz -107.37045 rang z=4,5
initial szz -87.84855 rang z=5,6
initial szz -68.32665 rang z=6,7
initial szz -48.80475 rang z=7,8
initial szz -29.28285 rang z=8,9
initial szz -9.76095 rang z=9,10

initial pp 98.1 rang z=-9.9999999999989E-05,9.9999999999989E-05
initial pp 88.29 rang z=0.9999,1.0001
prop density 1.49 rang z=0,1
initial pp 78.48 rang z=1.9999,2.0001
prop density 1.49 rang z=1,2
initial pp 68.67 rang z=2.9999,3.0001
prop density 1.49 rang z=2,3
initial pp 58.86 rang z=3.9999,4.0001
prop density 1.49 rang z=3,4
initial pp 49.05 rang z=4.9999,5.0001
prop density 1.49 rang z=4,5
initial pp 39.24 rang z=5.9999,6.0001
prop density 1.49 rang z=5,6
initial pp 29.43 rang z=6.9999,7.0001
prop density 1.49 rang z=6,7
initial pp 19.62 rang z=7.9999,8.0001
prop density 1.49 rang z=7,8
initial pp 9.81 rang z=8.9999,9.0001
prop density 1.49 rang z=8,9
initial pp 0 rang z=9.9999,10.0001
prop density 1.49 rang z=9,10
fix pp 0.0 rang z = 10,10

range name base z = 0
range name left x = 0
range name right x = 1

fix x range=nrange left
fix x range=nrange right
fix y
fix x y z range=nrange base
solve force 1

free x y z
fix y range= nrange base not

def tienodes
 local i = 0;
   loop i (1,10)
    global masternode   = gp_id(gp_near(0.0,0.0,0 + i*1))
    global slavenode1   = gp_id(gp_near(0.0,1.0,0 + i*1))
    global slavenode2   = gp_id(gp_near(1.0,1.0,0 + i*1))
    global slavenode3   = gp_id(gp_near(1.0,0.0,0 + i*1))
    Command
    attach gp @slavenode1 gp @masternode
    attach gp @slavenode2 gp @masternode
    attach gp @slavenode3 gp @masternode
    endcommand
  end_loop
end
@tienodes

;------------HISTORIES-----------------------------------------------------------------------------------;
hist unbal
hist dytime

;------------DEFINE SHEAR STRAIN-------------------------------------------------------------------------;
def shearStrain_z1
  array ai_z1(6)
  local pnt = find_zone(1)
  local dum = z_fsi(pnt, ai_z1)
  shearStrain_z1 = ai_z1(6) * 2;
End
def shearStrain_z2
  array ai_z2(6)
  local pnt = find_zone(2)
  local dum = z_fsi(pnt, ai_z2)
  shearStrain_z2 = ai_z2(6) * 2;
End
def shearStrain_z3
  array ai_z3(6)
  local pnt = find_zone(3)
  local dum = z_fsi(pnt, ai_z3)
  shearStrain_z3 = ai_z3(6) * 2;
End
def shearStrain_z4
  array ai_z4(6)
  local pnt = find_zone(4)
  local dum = z_fsi(pnt, ai_z4)
  shearStrain_z4 = ai_z4(6) * 2;
End
def shearStrain_z5
  array ai_z5(6)
  local pnt = find_zone(5)
  local dum = z_fsi(pnt, ai_z5)
  shearStrain_z5 = ai_z5(6) * 2;
End
def shearStrain_z6
  array ai_z6(6)
  local pnt = find_zone(6)
  local dum = z_fsi(pnt, ai_z6)
  shearStrain_z6 = ai_z6(6) * 2;
End
def shearStrain_z7
  array ai_z7(6)
  local pnt = find_zone(7)
  local dum = z_fsi(pnt, ai_z7)
  shearStrain_z7 = ai_z7(6) * 2;
End
def shearStrain_z8
  array ai_z8(6)
  local pnt = find_zone(8)
  local dum = z_fsi(pnt, ai_z8)
  shearStrain_z8 = ai_z8(6) * 2;
End
def shearStrain_z9
  array ai_z9(6)
  local pnt = find_zone(9)
  local dum = z_fsi(pnt, ai_z9)
  shearStrain_z9 = ai_z9(6) * 2;
End
def shearStrain_z10
  array ai_z10(6)
  local pnt = find_zone(10)
  local dum = z_fsi(pnt, ai_z10)
  shearStrain_z10 = ai_z10(6) * 2;
End

;------------RECORD EXCESS PORE PRESSURE----------------------------------------------------------------;
def inigpp
  array pp0(11)
  loop m(1,11)
    pp0(m) = gp_pp(gp_near(0, 0, m - 1))
  end_loop
End
@inigpp

def epp1
  epp1 = gp_pp(gp_near(0, 0, 1)) - pp0(1)
End

def epp2
  epp2 = gp_pp(gp_near(0, 0, 2)) - pp0(2)
End

def epp3
  epp3 = gp_pp(gp_near(0, 0, 3)) - pp0(3)
End

def epp4
  epp4 = gp_pp(gp_near(0, 0, 4)) - pp0(4)
End

def epp5
  epp5 = gp_pp(gp_near(0, 0, 5)) - pp0(5)
End

def epp6
  epp6 = gp_pp(gp_near(0, 0, 6)) - pp0(6)
End

def epp7
  epp7 = gp_pp(gp_near(0, 0, 7)) - pp0(7)
End

def epp8
  epp8 = gp_pp(gp_near(0, 0, 8)) - pp0(8)
End

def epp9
  epp9 = gp_pp(gp_near(0, 0, 9)) - pp0(9)
End

def epp10
  epp10 = gp_pp(gp_near(0, 0, 10)) - pp0(10)
End

def epp11
  epp11 = gp_pp(gp_near(0, 0, 11)) - pp0(11)
End


;------------DEFINE CONFINEMENT-------------------------------------------------------------------------;
def eff_stress_m1
  global pzm = z_near(0.5, 0.5, 0.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m1 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m2
  global pzm = z_near(0.5, 0.5, 1.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m2 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m3
  global pzm = z_near(0.5, 0.5, 2.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m3 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m4
  global pzm = z_near(0.5, 0.5, 3.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m4 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m5
  global pzm = z_near(0.5, 0.5, 4.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m5 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m6
  global pzm = z_near(0.5, 0.5, 5.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m6 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m7
  global pzm = z_near(0.5, 0.5, 6.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m7 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m8
  global pzm = z_near(0.5, 0.5, 7.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m8 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m9
  global pzm = z_near(0.5, 0.5, 8.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m9 = -1 * (tmp + z_pp(pzm))
End
def eff_stress_m10
  global pzm = z_near(0.5, 0.5, 9.5)
  local tmp = (z_sxx(pzm) + z_syy(pzm) + z_szz(pzm))/3.0
  eff_stress_m10 = -1 * (tmp + z_pp(pzm))
End

;------------SHEAR STRESS--------------------------------------------------------------------------------;
hist add zone sxz 0.5,0.5,0.5
hist add zone sxz 0.5,0.5,1.5
hist add zone sxz 0.5,0.5,2.5
hist add zone sxz 0.5,0.5,3.5
hist add zone sxz 0.5,0.5,4.5
hist add zone sxz 0.5,0.5,5.5
hist add zone sxz 0.5,0.5,6.5
hist add zone sxz 0.5,0.5,7.5
hist add zone sxz 0.5,0.5,8.5
hist add zone sxz 0.5,0.5,9.5

;------------SHEAR STRAIN--------------------------------------------------------------------------------;
hist add fish @shearStrain_z1
hist add fish @shearStrain_z2
hist add fish @shearStrain_z3
hist add fish @shearStrain_z4
hist add fish @shearStrain_z5
hist add fish @shearStrain_z6
hist add fish @shearStrain_z7
hist add fish @shearStrain_z8
hist add fish @shearStrain_z9
hist add fish @shearStrain_z10

;------------ACCELERATION--------------------------------------------------------------------------------;
hist add gp xacc 0,0,0
hist add gp xacc 0,0,1
hist add gp xacc 0,0,2
hist add gp xacc 0,0,3
hist add gp xacc 0,0,4
hist add gp xacc 0,0,5
hist add gp xacc 0,0,6
hist add gp xacc 0,0,7
hist add gp xacc 0,0,8
hist add gp xacc 0,0,9
hist add gp xacc 0,0,10

;------------DISPLACEMENT--------------------------------------------------------------------------------;
hist add gp xdisp 0,0,0
hist add gp xdisp 0,0,1
hist add gp xdisp 0,0,2
hist add gp xdisp 0,0,3
hist add gp xdisp 0,0,4
hist add gp xdisp 0,0,5
hist add gp xdisp 0,0,6
hist add gp xdisp 0,0,7
hist add gp xdisp 0,0,8
hist add gp xdisp 0,0,9
hist add gp xdisp 0,0,10

;------------VERTICAL STRESS-----------------------------------------------------------------------------;
hist add zone szz 0.5,0.5,0.5
hist add zone szz 0.5,0.5,1.5
hist add zone szz 0.5,0.5,2.5
hist add zone szz 0.5,0.5,3.5
hist add zone szz 0.5,0.5,4.5
hist add zone szz 0.5,0.5,5.5
hist add zone szz 0.5,0.5,6.5
hist add zone szz 0.5,0.5,7.5
hist add zone szz 0.5,0.5,8.5
hist add zone szz 0.5,0.5,9.5

;------------LATERAL STRESS------------------------------------------------------------------------------;
hist add zone sxx 0.5,0.5,0.5
hist add zone sxx 0.5,0.5,1.5
hist add zone sxx 0.5,0.5,2.5
hist add zone sxx 0.5,0.5,3.5
hist add zone sxx 0.5,0.5,4.5
hist add zone sxx 0.5,0.5,5.5
hist add zone sxx 0.5,0.5,6.5
hist add zone sxx 0.5,0.5,7.5
hist add zone sxx 0.5,0.5,8.5
hist add zone sxx 0.5,0.5,9.5

;------------PORE PRESSURE-------------------------------------------------------------------------------;
hist add gp pp 0,0,0
hist add gp pp 0,0,1
hist add gp pp 0,0,2
hist add gp pp 0,0,3
hist add gp pp 0,0,4
hist add gp pp 0,0,5
hist add gp pp 0,0,6
hist add gp pp 0,0,7
hist add gp pp 0,0,8
hist add gp pp 0,0,9
hist add gp pp 0,0,10

;------------CONFINEMENT-------------------------------------------------------------------------------;
hist add fish @eff_stress_m1
hist add fish @eff_stress_m2
hist add fish @eff_stress_m3
hist add fish @eff_stress_m4
hist add fish @eff_stress_m5
hist add fish @eff_stress_m6
hist add fish @eff_stress_m7
hist add fish @eff_stress_m8
hist add fish @eff_stress_m9
hist add fish @eff_stress_m10

;------------APPLICATION OF ACCELERATION-----------------------------------------------------------------;
set dyn on
set dyn damp=rayleigh 0.017322 2.7865
table 1 read SineWave.txt
apply xacc 2 his=table 1 yacc=0 zacc=0 range=nrange base

;------------DYNAMIC ANALYSIS----------------------------------------------------------------------------;
solve age 5

history write 3 vs 13 table StressStrain_at_Base
history write 8 vs 18 table StressStrain_at_Midheight
history write 12 vs 22 table StressStrain_at_Ground_Surface
history write 23 vs 2 table Acc_at_Base
history write 28 vs 2 table Acc_at_Midheight
history write 33 vs 2 table Acc_at_Ground_Surface
history write 34 vs 2 table Disp_at_Base
history write 39 vs 2 table Disp_at_Midheight
history write 44 vs 2 table Disp_at_Ground_Surface
history write 3 vs 76 table Conf_at_Base
history write 8 vs 81 table Conf_at_Midheight
history write 12 vs 85 table Conf_at_Ground_Surface
history write 65 vs 2 table PP_at_Base
history write 71 vs 2 table PP_at_Midheight
history write 75 vs 2 table PP_at_Ground_Surface

def substract_table(tab,tab2)
   local tbl_pnt = get_table(tab)
   local tbl_pnt2 = get_table(tab2)
   local tbl_size = table_size(tbl_pnt)
   local i = 0;
   loop i (1,tbl_size)
     local s = ytable(tbl_pnt,i)
     local t = ytable(tbl_pnt2,i)
     ytable(tbl_pnt, i) = s - t
   end_loop
End

def scalex_table(tab,scale)
   local tbl_pnt = get_table(tab)
   local tbl_size = table_size(tbl_pnt)
   local i = 0;
   loop i (1,tbl_size)
     local s = xtable(tbl_pnt,i)
     xtable(tbl_pnt, i) = s * scale
   end_loop
End

@scalex_table(2,100)
@scalex_table(3,100)
@scalex_table(4,100)
@substract_table(9,8)
@substract_table(10,8)

save Example
return
