;************************************************************************************     
;*                             					                                    *
;*                                                                                  *
;*                            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 'Qmk-Lama_del_Pacifico' slot 'terreno' 

;;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 1631.54 young 2.e6 poisson 0.35 range group  'Qmk-Lama_del_Pacifico' slot 'terreno'


;;DEFINICION DEL NIVEL FREATICO
zone water density 1000 ; kg
;;model gravity 10

;_____________________Water LEVEL__________________________________

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 'Qmk-Lama_del_Pacifico' 

;_________________________Initialisation Lateral Stresses Relation_____________________________
zone initialize-stresses ratio 0.6 total range group  'Qmk-Lama_del_Pacifico' slot 'terreno'

;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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'
;_______________________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_onlyLama.sav'