! ============================================================
! M09 - Fase de extraccion y validacion
! Hereda criterios de M07-M08. No escribe CSV (lo hace el driver).
! ============================================================

/POST1
*GET,n_sets,ACTIVE,0,SET,NSET
solution_available=0
*IF,n_sets,GE,1,THEN
  SET,LAST
  solution_available=1
*ENDIF
*IF,solution_available,NE,1,THEN
  /COM,ERROR: no result set available - inspect .err/.out before extracting
  /EOF
*ENDIF
RSYS,0

! --- Desplazamiento medio de la cara de punta ---
uy_sum=0
uy_tip_min=0
uy_tip_max=0
CMSEL,S,tip_nodes
node_id=0
*DO,j,1,n_tip
  node_id=NDNEXT(node_id)
  *GET,uy_node,NODE,node_id,U,Y
  uy_sum=uy_sum+uy_node
  *IF,j,EQ,1,THEN
    uy_tip_min=uy_node
    uy_tip_max=uy_node
  *ELSE
    *IF,uy_node,LT,uy_tip_min,THEN
      uy_tip_min=uy_node
    *ENDIF
    *IF,uy_node,GT,uy_tip_max,THEN
      uy_tip_max=uy_node
    *ENDIF
  *ENDIF
*ENDDO
uy_tip_avg=uy_sum/n_tip
uy_spread=uy_tip_max-uy_tip_min
uy_spread_ratio=ABS(uy_spread)/ABS(uy_tip_avg)
ALLSEL,ALL

! --- Referencia de Euler-Bernoulli ---
inertia=beam_b*beam_h**3/12
uy_ref=tip_force*beam_l**3/(3*young*inertia)
uy_error=ABS(uy_tip_avg-uy_ref)/ABS(uy_ref)

! --- Tension normal en seccion interior (x alineado a malla) ---
div_x_probe=NINT(beam_l/mesh_h)
x_index=NINT(0.2*div_x_probe)
*IF,x_index,LT,1,THEN
  x_index=1
*ENDIF
*IF,x_index,GT,div_x_probe,THEN
  x_index=div_x_probe
*ENDIF
x_section=x_index*beam_l/div_x_probe
moment_section=ABS(tip_force)*(beam_l-x_section)
sigma_ref=moment_section*(beam_h/2)/inertia

SELTOL,select_tol
NSEL,S,LOC,X,x_section
NSEL,R,LOC,Y,beam_h
CM,stress_top_nodes,NODE
*GET,n_top_section,NODE,0,COUNT

sx_top_sum=0
node_id=0
*DO,j,1,n_top_section
  node_id=NDNEXT(node_id)
  *GET,sx_node,NODE,node_id,S,X
  sx_top_sum=sx_top_sum+sx_node
*ENDDO
sx_top_avg=sx_top_sum/n_top_section
ALLSEL,ALL

NSEL,S,LOC,X,x_section
NSEL,R,LOC,Y,0
CM,stress_bottom_nodes,NODE
*GET,n_bottom_section,NODE,0,COUNT

sx_bottom_sum=0
node_id=0
*DO,j,1,n_bottom_section
  node_id=NDNEXT(node_id)
  *GET,sx_node,NODE,node_id,S,X
  sx_bottom_sum=sx_bottom_sum+sx_node
*ENDDO
sx_bottom_avg=sx_bottom_sum/n_bottom_section
ALLSEL,ALL
SELTOL,

sigma_fea=(ABS(sx_top_avg)+ABS(sx_bottom_avg))/2
stress_error=ABS(sigma_fea-sigma_ref)/sigma_ref
stress_symmetry=ABS(ABS(sx_top_avg)-ABS(sx_bottom_avg))/sigma_ref
stress_sign_pass=1
*IF,sx_top_avg,LE,0,THEN
  stress_sign_pass=0
*ENDIF
*IF,sx_bottom_avg,GE,0,THEN
  stress_sign_pass=0
*ENDIF

! --- Equilibrio ---
x_ref=0
y_ref=beam_h/2
z_ref=beam_b/2
rfx=0
rfy=0
rfz=0
rmx=0
rmy=0
rmz=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
  rfx=rfx+rfx_node
  rfy=rfy+rfy_node
  rfz=rfz+rfz_node
  rmx=rmx+ry*rfz_node-rz*rfy_node
  rmy=rmy+rz*rfx_node-rx*rfz_node
  rmz=rmz+rx*rfy_node-ry*rfx_node
*ENDDO
ALLSEL,ALL

external_mz=(beam_l-x_ref)*tip_force
force_error=ABS(rfy+tip_force)/ABS(tip_force)
moment_error=ABS(rmz+external_mz)/ABS(external_mz)

! --- Topologia esperada (misma regla que build) ---
div_x=NINT(beam_l/mesh_h)
div_y=NINT(beam_h/mesh_h)
div_z=NINT(beam_b/mesh_h)
*IF,div_x,LT,1,THEN
  div_x=1
*ENDIF
*IF,div_y,LT,1,THEN
  div_y=1
*ENDIF
*IF,div_z,LT,1,THEN
  div_z=1
*ENDIF
expected_nodes=(div_x+1)*(div_y+1)*(div_z+1)
expected_elements=div_x*div_y*div_z

! --- Veredicto de validacion (M07-M08) ---
validation_pass=1
*IF,n_nodes,NE,expected_nodes,THEN
  validation_pass=0
*ENDIF
*IF,n_elements,NE,expected_elements,THEN
  validation_pass=0
*ENDIF
*IF,n_tip,LE,0,THEN
  validation_pass=0
*ENDIF
*IF,solution_available,NE,1,THEN
  validation_pass=0
*ENDIF
*IF,n_sets,NE,1,THEN
  validation_pass=0
*ENDIF
*IF,uy_error,GE,displacement_tol,THEN
  validation_pass=0
*ENDIF
*IF,uy_spread_ratio,GE,spread_tol,THEN
  validation_pass=0
*ENDIF
*IF,n_top_section,LE,0,THEN
  validation_pass=0
*ENDIF
*IF,n_bottom_section,LE,0,THEN
  validation_pass=0
*ENDIF
*IF,stress_sign_pass,NE,1,THEN
  validation_pass=0
*ENDIF
*IF,stress_error,GE,stress_tol,THEN
  validation_pass=0
*ENDIF
*IF,stress_symmetry,GE,symmetry_tol,THEN
  validation_pass=0
*ENDIF
*IF,force_error,GE,equilibrium_tol,THEN
  validation_pass=0
*ENDIF
*IF,moment_error,GE,equilibrium_tol,THEN
  validation_pass=0
*ENDIF

! --- Veredicto de admisibilidad de diseno ---
feasible=0
*IF,validation_pass,EQ,1,THEN
  *IF,ABS(uy_tip_avg),LE,uy_limit,AND,sigma_fea,LE,stress_limit,THEN
    feasible=1
  *ENDIF
*ENDIF

/COM,M09 EXTRACT case=%case_id%: uy_error=%uy_error% feasible=%feasible% validation_pass=%validation_pass%