! ============================================================
! M01 - Sintaxis y parametros
! Proyecto conductor: viga maciza en voladizo
! Unidades coherentes: m, kg, s, N, Pa
! ============================================================

/CLEAR,START
/FILNAME,m01_parametric_beam,1
/TITLE,M01 - Viga parametrizada
/UNITS,SI

! --- 1. CONTRATO DE ENTRADA ---
beam_l=1.0
beam_h=0.10
beam_b=0.05
young=210E9
nu=0.30
tip_force=-1000
mesh_h=0.025
uy_limit=0.00025
stress_limit=150E6

! --- 2. VALIDACION TEMPRANA ---
*IF,beam_l,LE,0,THEN
  /COM,ERROR: beam_l must be positive
  /EOF
*ENDIF
*IF,beam_h,LE,0,THEN
  /COM,ERROR: beam_h must be positive
  /EOF
*ENDIF
*IF,beam_b,LE,0,THEN
  /COM,ERROR: beam_b must be positive
  /EOF
*ENDIF
*IF,young,LE,0,THEN
  /COM,ERROR: young must be positive
  /EOF
*ENDIF
*IF,mesh_h,LE,0,THEN
  /COM,ERROR: mesh_h must be positive
  /EOF
*ENDIF
*IF,mesh_h,GT,beam_h/2,THEN
  /COM,WARNING: mesh_h is larger than beam_h/2
*ENDIF

! --- 3. MAGNITUDES DERIVADAS Y REFERENCIAS ---
inertia=beam_b*beam_h**3/12
uy_ref=tip_force*beam_l**3/(3*young*inertia)
uy_ref_abs=ABS(uy_ref)
sigma_ref=ABS(tip_force)*beam_l*(beam_h/2)/inertia

! --- 4. PREDIMENSIONADO ANALITICO ---
*DIM,heights,ARRAY,4
heights(1)=0.08
heights(2)=0.10
heights(3)=0.12
heights(4)=0.14
selected_h=0

*CFOPEN,m01_design_table,csv
*VWRITE
('height_m,uy_ref_m,sigma_ref_pa,passes')
*DO,i,1,4
  trial_h=heights(i)
  trial_i=beam_b*trial_h**3/12
  trial_uy=ABS(tip_force)*beam_l**3/(3*young*trial_i)
  trial_stress=ABS(tip_force)*beam_l*(trial_h/2)/trial_i
  passes=0
  *IF,trial_uy,LE,uy_limit,THEN
    *IF,trial_stress,LE,stress_limit,THEN
      passes=1
      *IF,selected_h,EQ,0,THEN
        selected_h=trial_h
      *ENDIF
    *ENDIF
  *ENDIF
  *VWRITE,trial_h,trial_uy,trial_stress,passes
  (E16.8,',',E16.8,',',E16.8,',',F2.0)
*ENDDO
*CFCLOS

! --- 5. MODELO FEA DEL CASO BASE ---
/PREP7
ET,1,SOLID185
MP,EX,1,young
MP,PRXY,1,nu

BLOCK,0,beam_l,0,beam_h,0,beam_b
ESIZE,mesh_h
VMESH,ALL

NSEL,S,LOC,X,0
D,ALL,ALL,0
ALLSEL,ALL

NSEL,S,LOC,X,beam_l
*GET,n_tip,NODE,0,COUNT
*IF,n_tip,LE,0,THEN
  /COM,ERROR: no nodes found at beam tip
  /EOF
*ENDIF
F,ALL,FY,tip_force/n_tip
ALLSEL,ALL
FINISH

/SOLU
ANTYPE,STATIC
SOLVE
FINISH

/POST1
SET,LAST
PLNSOL,U,SUM
FINISH

! Valores esperados del caso base:
! inertia   = 4.1666667E-6 m^4
! uy_ref    = -3.8095238E-4 m
! sigma_ref = 1.2000000E7 Pa
! selected_h = 0.12 m