;************************************************************************************     
;*                             					                                    *
;*                                                                                  *
;*                            TEST MODEL                                            *
;*                         CARLOS MORA  (AECOM)                                     *
;*			         	      	                                                    *
;*		Effective vs Total stresses_Assessing settlements with tunnel face pressure *									                            *
;************************************************************************************
;
model new

zone import 'malla'

fish automatic-create on
customer-title-1 'AECOM'
customer-title-2 'www.aecom.com'

model large-strain off

zone group 'Fha-Rellenos_cohesivos' slot 'terreno' range pos-z 16.278 19.721 ;UPDATED
zone group 'Qmk-Lama_del_Pacifico' slot 'terreno' range pos-z 2.954 16.278
zone group 'Ov-2-Suelo_residual' slot 'terreno' range pos-z 0.170 2.954
zone group 'TI-3-w_Limolitas_calcareas_laminadas_Fm_La_Boca' slot 'terreno' range pos-z -5.788 0.170
zone group 'TI-2-s_Lutitas_limolitas_Fm_La_Boca' slot 'terreno' range pos-z -60 -5.788
;;RUTINA QUE MARCA LOS NODOS DE AVANCE ;NEW
zone gridpoint group 'avance' range group 'nucleo' slot '2'

;;MODELO CONSTITUTIVO Y PROPIEDADES DE CONTORNO.

;______________Boundary condition faces_______________________________
;
zone face skin; assign groups for boundary faces
zone gridpoint fix velocity-x range group 'East' or 'West'
zone gridpoint fix velocity-y range group 'North' or 'South'
zone gridpoint fix velocity-z range group 'Bottom'
zone gridpoint fix velocity-x range group 'Bottom'
zone gridpoint fix velocity-y range group 'Bottom'

model gravity 0,0,-9.81

;______________Ground properties______________________________________

zone cmodel assign elas ;range position-z @z_min1 @z_max

;______________________________NEW PROPERTY__________________________________________________
;Weak soil
zone property density 1834.86 young 4.e6 poisson 0.35 range group  'Fha-Rellenos_cohesivos' slot 'terreno'
zone property density 1631.54 young 2.e6 poisson 0.35 range group  'Qmk-Lama_del_Pacifico' slot 'terreno'
zone property density 1988.44 young 35.e6 poisson 0.35 range group  'Ov-2-Suelo_residual' slot 'terreno'

;Rock formation

zone property density 2243.37 young 32.e6 poisson 0.3 range group 'TI-3-w_Limolitas_calcareas_laminadas_Fm_La_Boca' slot 'terreno'
zone property density 2202.58 young 195.e6 poisson 0.3 range group 'TI-2-s_Lutitas_limolitas_Fm_La_Boca' slot 'terreno'

;;DEFINICION DEL NIVEL FREATICO
zone water density 1000 ; kg
zone water plane nor 0. 0. 1. o 0 0 19.721

model solve ratio 1e-6

;__________Assign materials p-h (Plastick hardening and h-b Hoek Brown)
zone cmodel assign p-h range group 'Fha-Rellenos_cohesivos' or 'Qmk-Lama_del_Pacifico' or 'Ov-2-Suelo_residual' or 'TI-3-w_Limolitas_calcareas_laminadas_Fm_La_Boca' or 'TI-2-s_Lutitas_limolitas_Fm_La_Boca'

;_________________________Initialisation Lateral Stresses Relation_____________________________
zone initialize-stresses ratio 0.5 total range group  'Fha-Rellenos_cohesivos' slot 'terreno'
zone initialize-stresses ratio 0.6 total range group  'Qmk-Lama_del_Pacifico' slot 'terreno'
zone initialize-stresses ratio 0.9 total range group  'Ov-2-Suelo_residual' slot 'terreno'
zone initialize-stresses ratio 1. total range group  'TI-3-w_Limolitas_calcareas_laminadas_Fm_La_Boca' slot 'terreno'
zone initialize-stresses ratio 1. total range group  'TI-2-s_Lutitas_limolitas_Fm_La_Boca' slot 'terreno'

;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
zone property density 1834.86 pressure-reference=100.e3 exponent=0.01 over-consolidation-ratio 1. flag-smallstrain on ...
 friction 27 cohesion 3.e3 stiffness-50-ref 4.e6 stiffness-ur-reference 12.e6 poisson=0.35 stiffness-0-reference 95.e6 strain-70 1.9e-4 range group 'Fha-Rellenos_cohesivos';
zone property density 1631.546 pressure-reference=100.e3 exponent=0.01 over-consolidation-ratio 1. flag-smallstrain on ...
  friction 23 cohesion 0.e3 stiffness-50-ref 2.e6 stiffness-ur-reference 6.e6 poisson=0.35 stiffness-0-reference 40.e6 strain-70 3.8e-4 range group 'Qmk-Lama_del_Pacifico'
zone property density 1988.447 pressure-reference=100.e3 exponent=0.01 over-consolidation-ratio 1. flag-smallstrain on ...
 friction 29 cohesion 18.e3 stiffness-50-ref 35.e6 stiffness-ur-reference 105.e6 poisson=0.33 stiffness-0-reference 220.0e6 strain-70 1.9e-4 range group 'Ov-2-Suelo_residual'

zone property density 2243.376 pressure-reference=100.e3 exponent=0.01 over-consolidation-ratio 1. flag-smallstrain off ...
 friction 31.5 cohesion 51.e3 stiffness-50-ref 32.e6 stiffness-ur-reference 96.e6 poisson=0.3 range group 'TI-3-w_Limolitas_calcareas_laminadas_Fm_La_Boca' ;stiffness-0-reference 900.0e6
zone property density 2202.587 pressure-reference=100.e3 exponent=0.01 over-consolidation-ratio 1. flag-smallstrain off ...
 friction 28.3 cohesion 57.e3 stiffness-50-ref 195.e6 stiffness-ur-reference 585.e6 poisson=0.3 range group 'TI-2-s_Lutitas_limolitas_Fm_La_Boca' ;stiffness-0-reference 3470.0e6
;_______________________Initial Condition_________________
model solve ratio 1e-6
model sav 'ini'
model res 'ini'

zone gridpoint initialize displacement-x 0
zone gridpoint initialize displacement-y 0
zone gridpoint initialize displacement-z 0
zone gridpoint initialize displacement-small-x 0
zone gridpoint initialize displacement-small-y 0
zone gridpoint initialize displacement-small-z 0
zone initialize state 0
;
model largestrain on
;;
;;**************** EXCAVACION SECUENCIAL *********************

fish define parm
rad_exc = 13.51/2  ;Excavation radius
pase = 2.  ;; NEW:Changing 2 to 1

x_min = 0.
x_max = 60.
y_min = 0.
y_max = 40.
;_______________________________________________________________
x_centt = 0. ; center of the excavation X
z_centt = 0. ; center of the excavation Y
x_centexc = x_min
y_cenexc = 0. 

ny1 = 40 ;; NEW:Number of zones for the tunnel excavation 
nr = 20 ;; Changed from 10 to 20, more dense mesh
end

@parm
zone gridpoint group 'avance' range group 'nucleo' slot '2'

fish define marca_cuts
  pnt = gp.head
  loop while pnt#null
   if gp.group(pnt,'default')='avance'
    yyp = gp.pos.y(pnt)
      gp.extra(pnt,1) = (yyp - y_min)/pase
   end_if
    pnt = gp.next(pnt)
  end_loop
end
@marca_cuts

fish define coloca_frente
  global cut
  global yya
  pnt = gp.head
  loop while pnt#null
    cutnd = gp.extra(pnt,1)
    if math.abs(cutnd - cut) < 1e-4 then
       gp.pos.y(pnt) = yya
    end_if
    pnt = gp.next(pnt)
  end_loop
end


fish define contorno_suelo 
  array persuelo(500,1000,4)  
  deltang = math.pi/float(2*nr)
  loop i(1,ny1 + 1)
    yy0 = y_min + (i - 1)
    loop j(1,2*2*nr);2*nr+1)
      angulo = math.pi/2. - (j - 1) * deltang
      xx0 = x_centexc + rad_exc * math.cos(angulo)
      zz0 = y_cenexc + rad_exc*math.sin(angulo)
      pnt = gp.near(xx0,yy0,zz0)
      persuelo(i,j,1) = pnt
    end_loop
  end_loop
end
@contorno_suelo

[nr = 20]
fish define perfila_excavacion ;NEW
    global i_frente
    i_frente2=i_frente-1
    inc_ang = math.pi/float(2*nr)
    loop j(1,2*2*nr);2*nr+1)
      angulo = math.pi/2. - (j - 1)*inc_ang
      pnt=persuelo(i_frente,j,1)
      pnt2=persuelo(i_frente2,j,1)
      gp.pos.x(pnt)=x_centt+(rad_exc+1e-4)*math.cos(angulo)
      gp.pos.z(pnt)=z_centt+(rad_exc+1e-4)*math.sin(angulo)
      gp.pos.x(pnt2)=x_centt+(rad_exc+1e-4)*math.cos(angulo)
      gp.pos.z(pnt2)=z_centt+(rad_exc+1e-4)*math.sin(angulo)
    end_loop
end

zone face group 'excavacion' slot 'lateral_pressure' internal range cyl e-1 0,-1,0 e-2 0,160,0 r 6.8 group 'roca'

fish define avanza_tunel
 p_fre=-231.77e3 ;Pressure at tunnel axis   

  loop cut(1,20)
    yyb=(y_min)+(cut-1)*pase
    yya=(y_min)+cut*pase
    yyam=yya-.01
    yyap=yya+.01
    i_frente=pase*cut+1 ;NEW positioning at the front face, this change is because there are double of zones in the aligment direction        ;Indice para entrar en los punteros del frente que estamos excavando
    command
     ;excavacion
    zone delete range group 'nucleo' slot '2' position-y @yyb @yya
;
    @coloca_frente
    @perfila_excavacion
    ;@contorno_permeable   ;
            
      zone face apply stress-normal @p_fre gradient 0. 0. 0.012e6 range group 'nucleo' slot '2' position-y @yyam @yyap
      zone face apply stress-normal @p_fre gradient 0. [20e3/14] 0.012e6 origin 0,[yya],0 range group 'excavacion' slot 'lateral_pressure' pos-y 0. [yya]

;_____________________Launching Calculation Stepping_______________
    model solve ratio 1e-6
    end_command

 end_loop
end

@avanza_tunel

model save 'EffectiveStresses_allmat.sav'