! ============================================================
! M10 - No linealidad geometrica y material
! Tres casos sobre viga esbelta: lineal, NLGEOM elastico, elastoplastico.
! Unidades coherentes: m, kg, s, N, Pa.
! ============================================================

/CLEAR,START
/FILNAME,m10_nonlinear,1
/TITLE,M10 - No linealidad geometrica y material
/UNITS,SI

beam_l=0.5
beam_h=0.04
beam_b=0.05
mesh_h=0.02
young=210E9
nu=0.30
tip_force=-10000
n_substeps=10
yield_stress=250E6
select_tol=MIN(beam_h,beam_b)*1E-5
equilibrium_tol=0.005
min_geom_dev=0.0005
max_geom_dev=0.25
min_plastic_dev=0.05
uy_tol=0.08

inertia=beam_b*beam_h**3/12
uy_ref=tip_force*beam_l**3/(3*young*inertia)

/PREP7                              ! CASO 1 - BASELINE LINEAL (NLGEOM OFF)
ET,1,SOLID185
KEYOPT,1,2,3
MP,EX,1,young
MP,PRXY,1,nu
TYPE,1
MAT,1
BLOCK,0,beam_l,0,beam_h,0,beam_b

div_x=beam_l/mesh_h
div_y=beam_h/mesh_h
div_z=beam_b/mesh_h
MSHAPE,0,3D
MSHKEY,1
LSEL,S,LENGTH,,beam_l
LESIZE,ALL,,,div_x,,1
LSEL,S,LENGTH,,beam_h
LESIZE,ALL,,,div_y,,1
LSEL,S,LENGTH,,beam_b
LESIZE,ALL,,,div_z,,1
ALLSEL,ALL
VMESH,ALL

*GET,n_nodes,NODE,0,COUNT
*GET,n_elements,ELEM,0,COUNT

SELTOL,select_tol
NSEL,S,LOC,X,0
CM,fixed_nodes,NODE
*GET,n_fixed,NODE,0,COUNT
NSEL,S,LOC,X,beam_l
CM,tip_nodes,NODE
*GET,n_tip,NODE,0,COUNT
NSEL,S,LOC,X,beam_l
NSEL,R,LOC,Y,0
NSEL,R,LOC,Z,0
*GET,corner_node,NODE,0,NUM,MIN
ALLSEL,ALL
SELTOL,

FINISH
/SOLU
CMSEL,S,fixed_nodes
D,ALL,ALL,0
ALLSEL,ALL
CMSEL,S,tip_nodes
F,ALL,FY,tip_force/n_tip
ALLSEL,ALL
ANTYPE,STATIC
NLGEOM,OFF                          ! NLGEOM,OFF: caso lineal de referencia sin grandes desplazamientos.
KBC,1
NSUBST,1                            ! NSUBST,1: un unico substep (carga instantanea, valido solo en lineal).
OUTRES,ALL,LAST                     ! OUTRES,ALL,LAST: guarda solo el ultimo set (suficiente con un substep).
SOLVE
FINISH

/POST1
SET,LAST
*GET,n_sets_linear,ACTIVE,0,SET,NSET
*GET,uy_linear,NODE,corner_node,U,Y

! Equilibrio caso lineal (dentro de POST1)
x_ref=0
y_ref=beam_h/2
z_ref=beam_b/2
rfy_linear=0
rmz_linear=0
CMSEL,S,fixed_nodes
node_id=0
*DO,j,1,n_fixed
  node_id=NDNEXT(node_id)
  *GET,node_x,NODE,node_id,LOC,X
  *GET,node_y,NODE,node_id,LOC,Y
  *GET,node_z,NODE,node_id,LOC,Z
  *GET,rfx_node,NODE,node_id,RF,FX
  *GET,rfy_node,NODE,node_id,RF,FY
  *GET,rfz_node,NODE,node_id,RF,FZ
  rx=node_x-x_ref
  ry=node_y-y_ref
  rz=node_z-z_ref
  rfy_linear=rfy_linear+rfy_node
  rmz_linear=rmz_linear+rx*rfy_node-ry*rfx_node
*ENDDO
ALLSEL,ALL
FINISH

external_mz=(beam_l-x_ref)*tip_force
force_error_linear=ABS(rfy_linear+tip_force)/ABS(tip_force)
moment_error_linear=ABS(rmz_linear+external_mz)/ABS(external_mz)
uy_error_linear=ABS(uy_linear-uy_ref)/ABS(uy_ref)

PARSAV,ALL,m10_state,par            ! PARSAV: guarda parametros escalares en archivo .par para el siguiente caso.

/CLEAR,NOSTART                      ! CASO 2 - NO LINEALIDAD GEOMETRICA (NLGEOM ON, ELASTICO)
/FILNAME,m10_geom,1
/TITLE,M10 - Caso geometrico NLGEOM ON
/UNITS,SI
PARRES,NEW,m10_state,par            ! PARRES: restaura parametros desde .par tras /CLEAR,NOSTART.
select_tol=MIN(beam_h,beam_b)*1E-5
div_x=beam_l/mesh_h
div_y=beam_h/mesh_h
div_z=beam_b/mesh_h

/PREP7
ET,1,SOLID185
KEYOPT,1,2,3
MP,EX,1,young
MP,PRXY,1,nu
TYPE,1
MAT,1
BLOCK,0,beam_l,0,beam_h,0,beam_b
LSEL,S,LENGTH,,beam_l
LESIZE,ALL,,,div_x,,1
LSEL,S,LENGTH,,beam_h
LESIZE,ALL,,,div_y,,1
LSEL,S,LENGTH,,beam_b
LESIZE,ALL,,,div_z,,1
ALLSEL,ALL
VMESH,ALL

SELTOL,select_tol
NSEL,S,LOC,X,0
CM,fixed_nodes,NODE
*GET,n_fixed2,NODE,0,COUNT
NSEL,S,LOC,X,beam_l
CM,tip_nodes,NODE
*GET,n_tip2,NODE,0,COUNT
NSEL,S,LOC,X,beam_l
NSEL,R,LOC,Y,0
NSEL,R,LOC,Z,0
*GET,corner_node2,NODE,0,NUM,MIN
ALLSEL,ALL
SELTOL,

FINISH
/SOLU
CMSEL,S,fixed_nodes
D,ALL,ALL,0
ALLSEL,ALL
CMSEL,S,tip_nodes
F,ALL,FY,tip_force/n_tip2
ALLSEL,ALL
ANTYPE,STATIC
NLGEOM,ON                           ! NLGEOM,ON: activa no linealidad geometrica (grandes desplazamientos).
AUTOTS,OFF
KBC,0
NSUBST,n_substeps,10,10             ! NSUBST: rampa de carga en n_substeps incrementos iguales.
NEQIT,50                            ! NEQIT: limite maximo de iteraciones de equilibrio por substep.
OUTRES,ALL,ALL                      ! OUTRES,ALL,ALL: guarda todos los substeps para trazar la curva de carga.
TIME,1
SOLVE
FINISH

/POST1
*GET,n_sets_geom,ACTIVE,0,SET,NSET
SET,LAST
*GET,uy_geom,NODE,corner_node2,U,Y

geom_dev_ratio=ABS((uy_geom-uy_linear)/uy_linear)

*CFOPEN,m10_load_curve,csv          ! Curva carga-desplazamiento (POST1)
*VWRITE
('substep,load_factor,uy_corner_m,case_tag')
*DO,i,1,n_sets_geom
  SET,1,1,,,i
  *GET,load_factor,ACTIVE,0,SET,TIME
  *GET,uy_step,NODE,corner_node2,U,Y
  case_tag=2
  *VWRITE,i,load_factor,uy_step,case_tag
  (F6.0,',',E16.8,',',E16.8,',',F3.0)
*ENDDO
*CFCLOS

! Equilibrio caso geometrico
rfy_geom=0
rmz_geom=0
SET,LAST
CMSEL,S,fixed_nodes
node_id=0
*DO,j,1,n_fixed2
  node_id=NDNEXT(node_id)
  *GET,node_x,NODE,node_id,LOC,X
  *GET,node_y,NODE,node_id,LOC,Y
  *GET,rfx_node,NODE,node_id,RF,FX
  *GET,rfy_node,NODE,node_id,RF,FY
  rx=node_x-x_ref
  ry=node_y-y_ref
  rfy_geom=rfy_geom+rfy_node
  rmz_geom=rmz_geom+rx*rfy_node-ry*rfx_node
*ENDDO
ALLSEL,ALL
FINISH

force_error_geom=ABS(rfy_geom+tip_force)/ABS(tip_force)
moment_error_geom=ABS(rmz_geom+external_mz)/ABS(external_mz)

PARSAV,ALL,m10_state,par            ! PARSAV: conserva estado numerico antes del caso elastoplastico.

/CLEAR,NOSTART                      ! CASO 3 - PLASTICIDAD BILINEAL (NLGEOM ON + TB,BISO)
/FILNAME,m10_plastic,1
/TITLE,M10 - Caso elastoplastico
/UNITS,SI
PARRES,NEW,m10_state,par            ! PARRES: recupera parametros del caso geometrico previo.
select_tol=MIN(beam_h,beam_b)*1E-5
div_x=beam_l/mesh_h
div_y=beam_h/mesh_h
div_z=beam_b/mesh_h

/PREP7
ET,1,SOLID185
KEYOPT,1,2,3
MP,EX,1,young
MP,PRXY,1,nu
TB,BISO,1,1                         ! TB,BISO: plasticidad bilineal isotropica (fluencia + endurecimiento).
TBDATA,1,yield_stress,0             ! TBDATA: limite de fluencia (yield_stress) y pendiente de endurecimiento (0=rigido-perfecto).
TYPE,1
MAT,1
BLOCK,0,beam_l,0,beam_h,0,beam_b
LSEL,S,LENGTH,,beam_l
LESIZE,ALL,,,div_x,,1
LSEL,S,LENGTH,,beam_h
LESIZE,ALL,,,div_y,,1
LSEL,S,LENGTH,,beam_b
LESIZE,ALL,,,div_z,,1
ALLSEL,ALL
VMESH,ALL

SELTOL,select_tol
NSEL,S,LOC,X,0
CM,fixed_nodes,NODE
*GET,n_fixed3,NODE,0,COUNT
NSEL,S,LOC,X,beam_l
CM,tip_nodes,NODE
*GET,n_tip3,NODE,0,COUNT
NSEL,S,LOC,X,beam_l
NSEL,R,LOC,Y,0
NSEL,R,LOC,Z,0
*GET,corner_node3,NODE,0,NUM,MIN
ALLSEL,ALL
SELTOL,

FINISH
/SOLU
CMSEL,S,fixed_nodes
D,ALL,ALL,0
ALLSEL,ALL
CMSEL,S,tip_nodes
F,ALL,FY,tip_force/n_tip3
ALLSEL,ALL
ANTYPE,STATIC
NLGEOM,ON                           ! NLGEOM,ON: mantiene grandes desplazamientos activos con plasticidad.
AUTOTS,OFF
KBC,0
NSUBST,n_substeps,10,10             ! NSUBST: misma rampa de carga que el caso geometrico elastico.
NEQIT,50                            ! NEQIT: tolera mas iteraciones cuando entra plasticidad en la iteracion.
OUTRES,ALL,ALL                      ! OUTRES,ALL,ALL: conserva todos los substeps para auditoria y curva.
TIME,1
SOLVE
FINISH

/POST1
*GET,n_sets_plastic,ACTIVE,0,SET,NSET
SET,LAST
*GET,uy_plastic,NODE,corner_node3,U,Y

plastic_dev_ratio=ABS((uy_plastic-uy_geom)/uy_geom)

! Equilibrio caso plastico
rfy_plastic=0
rmz_plastic=0
SET,LAST
CMSEL,S,fixed_nodes
node_id=0
*DO,j,1,n_fixed3
  node_id=NDNEXT(node_id)
  *GET,node_x,NODE,node_id,LOC,X
  *GET,node_y,NODE,node_id,LOC,Y
  *GET,rfx_node,NODE,node_id,RF,FX
  *GET,rfy_node,NODE,node_id,RF,FY
  rx=node_x-x_ref
  ry=node_y-y_ref
  rfy_plastic=rfy_plastic+rfy_node
  rmz_plastic=rmz_plastic+rx*rfy_node-ry*rfx_node
*ENDDO
ALLSEL,ALL
FINISH

force_error_plastic=ABS(rfy_plastic+tip_force)/ABS(tip_force)
moment_error_plastic=ABS(rmz_plastic+external_mz)/ABS(external_mz)

passes=1                            ! AUDITORIA Y CONTRATO
*IF,n_nodes,NE,312,THEN
  passes=0
*ENDIF
*IF,n_elements,NE,150,THEN
  passes=0
*ENDIF
*IF,n_sets_linear,NE,1,THEN
  passes=0
*ENDIF
*IF,n_sets_geom,NE,n_substeps,THEN
  passes=0
*ENDIF
*IF,n_sets_plastic,NE,n_substeps,THEN
  passes=0
*ENDIF
*IF,force_error_linear,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,force_error_geom,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,force_error_plastic,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,moment_error_linear,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,moment_error_geom,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,moment_error_plastic,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,uy_error_linear,GE,uy_tol,THEN
  passes=0
*ENDIF
*IF,geom_dev_ratio,LT,min_geom_dev,THEN
  passes=0
*ENDIF
*IF,geom_dev_ratio,GT,max_geom_dev,THEN
  passes=0
*ENDIF
*IF,plastic_dev_ratio,LT,min_plastic_dev,THEN
  passes=0
*ENDIF

*CFOPEN,m10_nonlinear_audit,csv
*VWRITE
('case,nlgeom,plasticity,n_substeps,n_nodes,n_elements,uy_corner_m,geom_dev_ratio,plastic_dev_ratio,force_error,moment_error,n_converged,passes')
*VWRITE,1,0,0,1,n_nodes,n_elements,uy_linear,0,0,force_error_linear,moment_error_linear,n_sets_linear,passes
(F3.0,',',F1.0,',',F1.0,',',F6.0,',',F10.0,',',F10.0,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',F6.0,',',F2.0)
*VWRITE,2,1,0,n_substeps,n_nodes,n_elements,uy_geom,geom_dev_ratio,0,force_error_geom,moment_error_geom,n_sets_geom,passes
(F3.0,',',F1.0,',',F1.0,',',F6.0,',',F10.0,',',F10.0,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',F6.0,',',F2.0)
*VWRITE,3,1,1,n_substeps,n_nodes,n_elements,uy_plastic,geom_dev_ratio,plastic_dev_ratio,force_error_plastic,moment_error_plastic,n_sets_plastic,passes
(F3.0,',',F1.0,',',F1.0,',',F6.0,',',F10.0,',',F10.0,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',F6.0,',',F2.0)
*CFCLOS

*CFOPEN,m10_summary,csv
*VWRITE
('uy_ref_m,uy_linear_m,uy_geom_m,uy_plastic_m,geom_dev_ratio,plastic_dev_ratio,study_passes')
*VWRITE,uy_ref,uy_linear,uy_geom,uy_plastic,geom_dev_ratio,plastic_dev_ratio,passes
(E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',F2.0)
*CFCLOS

*IF,passes,EQ,1,THEN
  /COM,M10 NONLINEAR STUDY PASSED
*ELSE
  /COM,M10 NONLINEAR STUDY FAILED - inspect m10_nonlinear_audit.csv and .out
  /STATUS,PARM
*ENDIF

FINISH