#Created by Zhaohui Yang 
#Modified by Jinchi Lu (jinlu@ucsd.edu)
#fully coupled, u-p formulation
#
#plane strain,  shear-beam type mesh with single material,  
#dynamic analysis,  English units (in, s, lbf)
#
#
wipe

# Unit system in this file: English

# define UNITS ---------------------------------------------
# Modified from Silvia Mazzoni & Frank McKenna, 2006
set in 1.;              # define basic units -- output units
set lb 1.;              # define basic units -- output units
set sec 1.;             # define basic units -- output units
set ft [expr 12.*$in];          # define engineering units
set ft2 [expr $ft*$ft];                 
set ft4 [expr $ft2*$ft2]; 
set kip [expr 1000.*$lb];
set ksi [expr $kip/pow($in,2)];
set psi [expr $ksi/1000.];
set lbf [expr $psi*$in*$in];                    # pounds force
set pcf [expr $lbf/pow($ft,3)];         # pounds per cubic foot
set psf [expr $lbf/pow($ft,2)];         # pounds per square foot
set in2 [expr $in*$in];                                 # inch^2
set in4 [expr $in*$in*$in*$in];                 # inch^4
set g [expr 32.2*$ft/pow($sec,2)];      # gravitational acceleration

set psfTokPa   4.788000e-002
set psiTokPa   6.894757e+000
set ksiTokPa   6.894757e+003
set ftTom   3.048000e-001
set inTom   2.540000e-002
set kipsTokN   4.448222e+000
set lbsTokN 4.448222e-3
set kipsftTokNm   1.355818e+000
set pcfToton_m3   1.601800e-002

set kPaTopsf  [expr 1./$psfTokPa];
set kPaTopsi  [expr 1./$psiTokPa];
set kPaToksi  [expr 1./$ksiTokPa];
set mToft [expr 1./$ftTom];
set mToin [expr 1./$inTom];
set kNTokips [expr 1./$kipsTokN];  
set kNTolbs [expr 1./$lbsTokN];
set kNmTokipsft  [expr 1./$kipsftTokNm];
set ton_m3Topcf  [expr 1./$pcfToton_m3];
set kkg_m3Toslug_ft3 1.94055

#
#some user defined variables
# 
set matOpt   1      ;# 1 = pressure depend; 
                    ;# 2 = pressure independ; 
set fmass [expr 1.*$kkg_m3Toslug_ft3/$ft4];
set smass [expr 2.*$kkg_m3Toslug_ft3/$ft4];
set G      [expr 6.e4*$kPaToksi*$ksi]     ;
set B     [expr 2.4e5*$kPaToksi*$ksi] ;
set bulk    [expr 2.2e6*$kPaToksi*$ksi] ;#fluid-solid combined bulk modulus
set vperm  [expr 5.e-4*$mToin]  ;#vertical permeability (m/s)
set hperm  [expr 5.e-4*$mToin]  ;#horizontal permeability (m/s)
set accGravity  $g    ;#acceleration of gravity
set vperm [expr $vperm/$accGravity/$fmass]  ;#actual value used in computation
set hperm [expr $hperm/$accGravity/$fmass]  ;#actual value used in computation
set press   [expr 0.*$kPaToksi*$ksi]    ;# isotropic consolidation pressure on quad element(s)
set loadBias 0.07     ;# Static shear load, in percentage of gravity load (=sin(inclination angle))

set accMul  [expr 2*$mToin];  ;# acc. multiplier
set accNam  whatever.acc  ;# file name for input acc. record 
set accDt   0.0166   ;# dt of input acc. record
set period   1.0      ;# Period for applied Sine wave
set deltaT   0.01    ;# time step for analysis
set numSteps 2400    ;# number of time steps
set gamma    0.6    ;# Newmark integration parameter

set massProportionalDamping   0. ;
set InitStiffnessProportionalDamping 0.002;

set numXele 1      ;# number of elements in x (H) direction
set numYele 10      ;# number of elements in y (V) direction
set xSize [expr 1.0*$mToin]       ;# x direction element size
set ySize [expr 1.0*$mToin]       ;# y direction element size

#############################################################
# BUILD MODEL

#create the ModelBuilder
model basic -ndm 2 -ndf 3

# define material and properties
switch $matOpt {
  1 {
    nDMaterial PressureDependMultiYield 1 2 $smass $G $B  31.4 .1 [expr 80*$kPaToksi*$ksi] 0.5 \
                                        26.5 0.1 .2 5 [expr 10*$kPaToksi*$ksi] 0.015 1.0 20 \
                                        0.6 0.9 0.02 0.7 [expr 101*$kPaToksi*$ksi] [expr 0.3*$kPaToksi*$ksi] 

  }
  2 {
    nDMaterial PressureIndependMultiYield 2 2 [expr 1.8*$kkg_m3Toslug_ft3/$ft4] [expr 4.e4*$kPaToksi*$ksi] [expr 2.e5*$kPaToksi*$ksi] [expr 40.*$kPaToksi*$ksi] .1 \
                                     0. [expr 100*$kPaToksi*$ksi] 0. 20
  }
}

set gravY -$accGravity                 ;#calc. gravity
set gravX [expr -$gravY*$loadBias]

# define nodes
set numXnode  [expr $numXele+1]
set numYnode  [expr $numYele+1]

for {set i 1} {$i <= $numXnode} {incr i 1} {
  for {set j 1} {$j <= $numYnode} {incr j 1} {
    set xdim [expr ($i-1)*$xSize]
    set ydim [expr ($j-1)*$ySize] 
    set nodeNum [expr $i + ($j-1)*$numXnode] 
    node $nodeNum $xdim $ydim 
  }
}

# define elements
for {set i 1} {$i <= $numXele} {incr i 1} {
  for {set j 1} {$j <= $numYele} {incr j 1} {
    set eleNum [expr $i + ($j-1)*$numXele] 
    set n1  [expr $i + ($j-1)*$numXnode] 
    set n2  [expr $i + ($j-1)*$numXnode + 1] 
    set n4  [expr $i + $j*$numXnode + 1] 
    set n3  [expr $i + $j*$numXnode] 
                                    #      thick maTag  bulk  mDensity  perm1  perm2  gravity      press   
    element quadUP $eleNum $n1 $n2 $n4 $n3 [expr 1.0*$mToin]  $matOpt $bulk $fmass  $hperm  $vperm $gravX $gravY $press   
  }
}  

#set material to elastic for gravity loading
updateMaterialStage -material $matOpt -stage 0  

# fix the base, and free surface drainage
for {set i 1} {$i <= $numXnode} {incr i 1} {
  fix $i 1 1 0
  set surfnode [expr ($numYnode-1)*$numXnode + $i] 
  fix $surfnode 0 0 1
}

# tie all disp. DOFs at same level
for {set i 1} {$i < $numYnode} {incr i 1} {
  set nodeNum1 [expr $i*$numXnode + 1]
  for {set j 2} {$j <= $numXnode} {incr j 1} {
    set nodeNum2 [expr $i*$numXnode + $j]
    equalDOF $nodeNum1 $nodeNum2 1 2 
  }
}

#############################################################
# GRAVITY APPLICATION (elastic behavior)

# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer
numberer RCM
system ProfileSPD
test NormDispIncr 1.0e-8 25 2
algorithm Newton
constraints Penalty 1.e18 1.e18
integrator Newmark 1.5  1.  
analysis Transient 


analyze 3 5.e3

# switch material stage from elastic (gravity) to plastic
switch $matOpt {
  1 {
     updateMaterialStage -material $matOpt -stage 1
  }
  2 {
     updateMaterialStage -material $matOpt -stage 1
  }
}
analyze 5 5.e3
# rezero time
wipeAnalysis
setTime 0.0
#loadConst -time 0.0
#############################################################
# NOW APPLY LOADING SEQUENCE AND ANALYZE (plastic)

#   base input motion                     
pattern UniformExcitation  1  1    -accel "Sine 0. 10. $period -factor $accMul"

#input motion through data file
#pattern UniformExcitation  1  1  -accel "Series -factor $accMul -filePath $accNam -dt $accDt"

#recorder for nodal variables along the vertical center line.
set nodeList {}
for {set i 0} {$i < $numYnode} {incr i 1} {
  lappend nodeList [expr $numXnode/2 + $i*$numXnode]
}

#define recorders for disp., excess pore pressure, and acc.
#Note: disp and acc outputs are relative to the base
eval "recorder Node -file disp  -time -node $nodeList -dof 1 2 -dT $deltaT disp"
eval "recorder Node -file pwp  -time -node $nodeList -dof 3 -dT $deltaT vel"
eval "recorder Node -file acc  -time -node $nodeList -dof 1 2 -dT $deltaT accel"

#stress/strain output at Gauss point 1 of each element along center line
set name1 "stress";   set name2 "strain";  
for {set i 1} {$i < $numYnode} {incr i 1} {
  set ele [expr $numXele-$numXele/2+($i-1)*$numXele] 
  set name11 [join [list $name1 $i] {}]
  set name21 [join [list $name2 $i] {}] 
  recorder Element -ele $ele  -time -file $name11  -dT $deltaT material 1 stress
  recorder Element -ele $ele  -time -file $name21  -dT $deltaT material 1 strain
}

#analysis options
constraints Penalty 1.e18 1.e18
test NormDispIncr 1.e-5 25 0
numberer RCM
algorithm Newton
system ProfileSPD
#some mass proportional and initial-stiffness proportional damping
rayleigh $massProportionalDamping 0.0 $InitStiffnessProportionalDamping 0.0
integrator Newmark $gamma  [expr pow($gamma+0.5, 2)/4]  
analysis VariableTransient 

#analyze 
set startT [clock seconds]
analyze $numSteps $deltaT [expr $deltaT/100] $deltaT 15
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."

wipe  #flush ouput stream