! ============================================================
! M07 - Postprocesado y validacion
! Visualiza, extrae y valida una viga SOLID185.
! Unidades coherentes: m, kg, s, N, Pa.
! ============================================================

/CLEAR,START
/FILNAME,m07_validation,1
/TITLE,M07 - Postprocesado y validacion
/UNITS,SI

case_id=0
beam_l=1.0
beam_h=0.10
beam_b=0.05
mesh_h=0.05
young=210E9
nu=0.30
tip_force=-1000
select_tol=MIN(beam_h,beam_b)*1E-5
equilibrium_tol=0.005
displacement_tol=0.05
spread_tol=0.01
trace_tol=0.001
stress_tol=0.10
symmetry_tol=0.05

! --- 1. MODELO, MALLA Y REGIONES HEREDADOS ---
/PREP7
ET,1,SOLID185
! KEYOPT(2)=3: simplified enhanced strain para flexion (recomendado por MAPDL).
! KEYOPT(2)=0 produce ~22 % de error en S,X con esta malla.
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
CSYS,0
NSEL,S,LOC,X,0
CM,fixed_nodes,NODE
*GET,n_fixed,NODE,0,COUNT
ALLSEL,ALL

NSEL,S,LOC,X,beam_l
CM,tip_nodes,NODE
*GET,n_tip,NODE,0,COUNT
ALLSEL,ALL
SELTOL,

! --- 2. CONDICIONES Y SOLUCION HEREDADAS DE M05-M06 ---
FINISH
/SOLU
CMSEL,S,fixed_nodes
D,ALL,ALL,0
ALLSEL,ALL

CMSEL,S,tip_nodes
force_per_node=tip_force/n_tip
F,ALL,FY,force_per_node
ALLSEL,ALL

ANTYPE,STATIC
NLGEOM,OFF
KBC,1
NSUBST,1
OUTRES,ALL,LAST
SOLVE
FINISH

! --- 3. ACTIVAR CONSCIENTEMENTE EL SET ---
/POST1
*GET,n_sets,ACTIVE,0,SET,NSET
solution_available=0
*IF,n_sets,GE,1,THEN
  SET,LAST
  solution_available=1
*ENDIF
RSYS,0

! --- 4. 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

! Nodo de esquina usado por M04: trazabilidad, no metrica principal.
SELTOL,select_tol
NSEL,S,LOC,X,beam_l
NSEL,R,LOC,Y,0
NSEL,R,LOC,Z,0
*GET,corner_node,NODE,0,NUM,MIN
*GET,uy_tip_corner,NODE,corner_node,U,Y
ALLSEL,ALL
SELTOL,

! Referencia de esquina con KEYOPT(2)=3 y mesh_h=0.05 (MAPDL Student 2025 R2).
uy_corner_ref=-3.79470E-4
corner_trace_error=ABS(uy_tip_corner-uy_corner_ref)/ABS(uy_corner_ref)

! --- 5. 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)

! --- 6. TENSION NORMAL EN UNA SECCION INTERIOR ---
x_section=0.2*beam_l
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

! --- 7. EQUILIBRIO COMO PUERTA DE ENTRADA ---
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)

! --- 8. CONTRATO DE ACEPTACION ---
passes=1
*IF,n_nodes,NE,126,THEN
  passes=0
*ENDIF
*IF,n_elements,NE,40,THEN
  passes=0
*ENDIF
*IF,n_tip,NE,6,THEN
  passes=0
*ENDIF
*IF,solution_available,NE,1,THEN
  passes=0
*ENDIF
*IF,n_sets,NE,1,THEN
  passes=0
*ENDIF
*IF,uy_error,GE,displacement_tol,THEN
  passes=0
*ENDIF
*IF,uy_spread_ratio,GE,spread_tol,THEN
  passes=0
*ENDIF
*IF,corner_trace_error,GE,trace_tol,THEN
  passes=0
*ENDIF
*IF,n_top_section,NE,2,THEN
  passes=0
*ENDIF
*IF,n_bottom_section,NE,2,THEN
  passes=0
*ENDIF
*IF,stress_sign_pass,NE,1,THEN
  passes=0
*ENDIF
*IF,stress_error,GE,stress_tol,THEN
  passes=0
*ENDIF
*IF,stress_symmetry,GE,symmetry_tol,THEN
  passes=0
*ENDIF
*IF,force_error,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,moment_error,GE,equilibrium_tol,THEN
  passes=0
*ENDIF

*CFOPEN,m07_validation_audit,csv
*VWRITE
('case,beam_l,mesh_h,n_nodes,n_elements,n_sets,n_tip,uy_tip_avg,uy_tip_corner,uy_ref,uy_error,sigma_fea,sigma_ref,stress_error,force_error,moment_error,passes')
*VWRITE,case_id,beam_l,mesh_h,n_nodes,n_elements,n_sets,n_tip,uy_tip_avg,uy_tip_corner,uy_ref,uy_error,sigma_fea,sigma_ref,stress_error,force_error,moment_error,passes
(F3.0,',',E16.8,',',E16.8,',',F10.0,',',F10.0,',',F6.0,',',F6.0,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',F2.0)
*CFCLOS

! --- 9. VISUALIZAR NO SUSTITUYE A EXTRAER ---
/GRAPHICS,POWER
PLDISP,2
PLNSOL,U,Y
PLNSOL,S,X
PLESOL,S,X

*IF,passes,EQ,1,THEN
  /COM,M07 VALIDATION PASSED
*ELSE
  /COM,M07 VALIDATION FAILED - inspect CSV and assumptions
*ENDIF

FINISH