! ============================================================
! M04 - Mallado, calidad y convergencia
! Familia de tres mallas hexaedricas estructuradas.
! D, F, SOLVE y /POST1 forman un banco de ensayo heredado.
! ============================================================

/CLEAR,START
/FILNAME,m04_mesh_study,1
/TITLE,M04 - Estudio de convergencia de malla
/UNITS,SI

beam_l=1.0
beam_h=0.10
beam_b=0.05
young=210E9
nu=0.30
tip_force=-1000
conv_tol=2.0
select_tol=MIN(beam_h,beam_b)*1E-5
n_cases=3

*DIM,mesh_size,ARRAY,n_cases
*DIM,div_x,ARRAY,n_cases
*DIM,div_y,ARRAY,n_cases
*DIM,div_z,ARRAY,n_cases
*DIM,n_nodes,ARRAY,n_cases
*DIM,n_elems,ARRAY,n_cases
*DIM,n_type1,ARRAY,n_cases
*DIM,n_mat1,ARRAY,n_cases
*DIM,n_shape_bad,ARRAY,n_cases
*DIM,mesh_pass,ARRAY,n_cases
*DIM,uy_tip,ARRAY,n_cases
*DIM,change_pct,ARRAY,n_cases
*DIM,selected,ARRAY,n_cases

mesh_size(1)=0.0500
mesh_size(2)=0.0250
mesh_size(3)=0.0125

! --- MODELO CONSTANTE PARA LOS TRES CASOS ---
/PREP7
ET,1,SOLID185
KEYOPT,1,2,0
MP,EX,1,young
MP,PRXY,1,nu
TYPE,1
MAT,1
BLOCK,0,beam_l,0,beam_h,0,beam_b

MSHAPE,0,3D
MSHKEY,1
SHPP,DEFAULT
SHPP,ON
SELTOL,select_tol

*DO,i,1,n_cases
  ! Las tres dimensiones son divisibles exactamente en el caso base.
  ! NX, NY y NZ son funciones intrinsecas de MAPDL:
  ! no pueden utilizarse como nombres de arrays.
  div_x(i)=beam_l/mesh_size(i)
  div_y(i)=beam_h/mesh_size(i)
  div_z(i)=beam_b/mesh_size(i)

  ! ESIZE fija el control global. LESIZE documenta y fuerza
  ! las divisiones en cada familia de lineas.
  ESIZE,mesh_size(i)

  LSEL,S,LENGTH,,beam_l
  LESIZE,ALL,,,div_x(i),,1
  LSEL,S,LENGTH,,beam_h
  LESIZE,ALL,,,div_y(i),,1
  LSEL,S,LENGTH,,beam_b
  LESIZE,ALL,,,div_z(i),,1
  ALLSEL,ALL

  VMESH,ALL
  ALLSEL,ALL

  ! --- PASAPORTE DE MALLA ---
  *GET,n_nodes(i),NODE,0,COUNT
  *GET,n_elems(i),ELEM,0,COUNT

  ESEL,S,TYPE,,1
  *GET,n_type1(i),ELEM,0,COUNT
  ALLSEL,ALL
  ESEL,S,MAT,,1
  *GET,n_mat1(i),ELEM,0,COUNT
  ALLSEL,ALL

  ! SUMMARY lista las metricas. CHECK conserva solo incidencias.
  SHPP,SUMMARY
  CHECK,ESEL,WARN
  *GET,n_shape_bad(i),ELEM,0,COUNT
  ALLSEL,ALL

  mesh_pass(i)=1
  expected_elems=div_x(i)*div_y(i)*div_z(i)
  expected_nodes=(div_x(i)+1)*(div_y(i)+1)*(div_z(i)+1)
  *IF,n_elems(i),NE,expected_elems,THEN
    mesh_pass(i)=0
  *ENDIF
  *IF,n_nodes(i),NE,expected_nodes,THEN
    mesh_pass(i)=0
  *ENDIF
  *IF,n_type1(i),NE,n_elems(i),THEN
    mesh_pass(i)=0
  *ENDIF
  *IF,n_mat1(i),NE,n_elems(i),THEN
    mesh_pass(i)=0
  *ENDIF
  *IF,n_shape_bad(i),NE,0,THEN
    mesh_pass(i)=0
  *ENDIF

  ! --- BANCO DE ENSAYO HEREDADO DE M00 ---
  ! M05-M07 explicaran formalmente cargas, solucion y resultados.
  NSEL,S,LOC,X,0
  CM,fixed_nodes,NODE
  D,ALL,ALL,0
  ALLSEL,ALL

  NSEL,S,LOC,X,beam_l
  *GET,n_tip,NODE,0,COUNT
  CM,tip_nodes,NODE
  force_per_node=tip_force/n_tip
  F,ALL,FY,force_per_node

  ! Nodo de medida definido geometricamente, nunca por ID.
  NSEL,R,LOC,Y,0
  NSEL,R,LOC,Z,0
  *GET,probe_node,NODE,0,NUM,MIN
  ALLSEL,ALL

  FINISH
  /SOLU
  ANTYPE,STATIC,NEW
  SOLVE
  FINISH

  /POST1
  SET,LAST
  *GET,uy_tip(i),NODE,probe_node,U,Y
  FINISH

  ! Conservar la ultima malla; limpiar entre los casos anteriores.
  *IF,i,LT,n_cases,THEN
    /PREP7
    ALLSEL,ALL
    DDELE,ALL,ALL
    FDELE,ALL,ALL
    CMDELE,fixed_nodes
    CMDELE,tip_nodes
    VCLEAR,ALL
  *ENDIF
*ENDDO

SELTOL,

! --- CONVERGENCIA: CADA CASO SE COMPARA CON EL SIGUIENTE ---
change_pct(1)=ABS(uy_tip(2)-uy_tip(1))/ABS(uy_tip(2))*100
change_pct(2)=ABS(uy_tip(3)-uy_tip(2))/ABS(uy_tip(3))*100
change_pct(3)=-1

selected_case=0
*IF,change_pct(1),LT,conv_tol,THEN
  *IF,mesh_pass(1),EQ,1,THEN
    *IF,mesh_pass(2),EQ,1,THEN
      selected_case=1
    *ENDIF
  *ENDIF
*ENDIF

*IF,selected_case,EQ,0,THEN
  *IF,change_pct(2),LT,conv_tol,THEN
    *IF,mesh_pass(2),EQ,1,THEN
      *IF,mesh_pass(3),EQ,1,THEN
        selected_case=2
      *ENDIF
    *ENDIF
  *ENDIF
*ENDIF

*IF,selected_case,GT,0,THEN
  selected(selected_case)=1
*ENDIF

*CFOPEN,m04_mesh_study,csv
*VWRITE
('case,mesh_h,nx,ny,nz,n_nodes,n_elements,n_type1,n_mat1,uy_tip,change_to_next_pct,mesh_pass,selected')
*DO,i,1,n_cases
  ! Copiar a escalares evita la vectorizacion automatica de *VWRITE.
  out_case=i
  out_h=mesh_size(i)
  out_dx=div_x(i)
  out_dy=div_y(i)
  out_dz=div_z(i)
  out_nodes=n_nodes(i)
  out_elems=n_elems(i)
  out_type=n_type1(i)
  out_mat=n_mat1(i)
  out_uy=uy_tip(i)
  out_change=change_pct(i)
  out_pass=mesh_pass(i)
  out_selected=selected(i)
  *VWRITE,out_case,out_h,out_dx,out_dy,out_dz,out_nodes,out_elems,out_type,out_mat,out_uy,out_change,out_pass,out_selected
  (F3.0,',',E16.8,',',F8.0,',',F8.0,',',F8.0,',',F12.0,',',F12.0,',',F12.0,',',F12.0,',',E16.8,',',E16.8,',',F2.0,',',F2.0)
*ENDDO
*CFCLOS

*IF,selected_case,GT,0,THEN
  /COM,M04 STUDY PASSED - first acceptable mesh selected
*ELSE
  /COM,M04 STUDY INCONCLUSIVE - add a finer mesh
*ENDIF

FINISH