! ============================================================
! M12 - Union atornillada idealizada
! Dos placas SOLID185 unidas por nodos piloto + CERIG + conector 1D.
! Unidades coherentes: m, kg, s, N, Pa.
! ============================================================
/CLEAR,START                        ! Casos: MPC184 (viga rigida) y BEAM188 (rigidez elastica).
/FILNAME,m12_bolt,1
/TITLE,M12 - Union atornillada idealizada
/UNITS,SI

plate_x=0.08
plate_z=0.04
plate_t=0.012
interface_gap=0.0001
mesh_h=0.004
bolt_x=0.04
bolt_z=0.02
head_absorb_r=0.010
thread_absorb_r=0.008
bolt_d=0.012
young=210E9
nu=0.30
fx_total=4000
fy_total=-3000
select_tol=mesh_h*1E-4
equilibrium_tol=0.005
force_match_tol=0.08
min_head_dep=4
min_thread_dep=4

/PREP7                              ! CASO 1 - MPC184 RIGID BEAM (KEYOPT 1=1, 2=1)
ET,1,SOLID185
KEYOPT,1,2,0
ET,2,MPC184                         ! MPC184: conector multipunto rigido entre nodos piloto de cabecera y rosca.
KEYOPT,2,1,1                        ! KEYOPT(1)=1: modo viga rigida (rigid beam).
KEYOPT,2,2,1                        ! KEYOPT(2)=1: eliminacion de DOF dependientes (expone SMISC en POST1).
MP,EX,1,young
MP,PRXY,1,nu
MP,EX,2,young
MP,PRXY,2,nu
TYPE,1
MAT,1
BLOCK,0,plate_x,0,plate_t,0,plate_z
BLOCK,0,plate_x,plate_t+interface_gap,2*plate_t+interface_gap,0,plate_z

div_x=MAX(1,NINT(plate_x/mesh_h))
div_y=MAX(1,NINT(plate_t/mesh_h))
div_z=MAX(1,NINT(plate_z/mesh_h))
MSHAPE,0,3D
MSHKEY,1
LSEL,S,LENGTH,,plate_x
LESIZE,ALL,,,div_x,,1
LSEL,S,LENGTH,,plate_t
LESIZE,ALL,,,div_y,,1
LSEL,S,LENGTH,,plate_z
LESIZE,ALL,,,div_z,,1
ALLSEL,ALL
VMESH,ALL

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

SELTOL,select_tol
NSEL,S,LOC,Y,0
CM,bottom_nodes,NODE
*GET,n_bottom,NODE,0,COUNT
NSEL,S,LOC,Y,2*plate_t+interface_gap
CM,top_load_nodes,NODE
*GET,n_top_load,NODE,0,COUNT
ALLSEL,ALL
SELTOL,

SELTOL,select_tol
NSEL,S,LOC,Y,2*plate_t+interface_gap
NSEL,R,LOC,X,bolt_x-head_absorb_r,bolt_x+head_absorb_r
NSEL,R,LOC,Z,bolt_z-head_absorb_r,bolt_z+head_absorb_r
CM,head_dep,NODE
*GET,n_head_dep,NODE,0,COUNT
CMSEL,S,head_dep
NSEL,R,LOC,X,bolt_x-select_tol,bolt_x+select_tol
NSEL,R,LOC,Z,bolt_z-select_tol,bolt_z+select_tol
*GET,head_pilot,NODE,0,NUM,MIN
CMSEL,S,head_dep
CERIG,head_pilot,ALL
ALLSEL,ALL

NSEL,S,LOC,Y,plate_t
NSEL,R,LOC,X,bolt_x-thread_absorb_r,bolt_x+thread_absorb_r
NSEL,R,LOC,Z,bolt_z-thread_absorb_r,bolt_z+thread_absorb_r
CM,thread_dep,NODE
*GET,n_thread_dep,NODE,0,COUNT
CMSEL,S,thread_dep
NSEL,R,LOC,X,bolt_x-select_tol,bolt_x+select_tol
NSEL,R,LOC,Z,bolt_z-select_tol,bolt_z+select_tol
*GET,thread_pilot,NODE,0,NUM,MIN
CMSEL,S,thread_dep
CERIG,thread_pilot,ALL
ALLSEL,ALL
SELTOL,

TYPE,2
REAL,1
MAT,2
E,head_pilot,thread_pilot
*GET,conn_elem,ELEM,0,NUM,MAX
ALLSEL,ALL
FINISH

/SOLU
ANTYPE,STATIC
NLGEOM,OFF
CMSEL,S,bottom_nodes
D,ALL,ALL,0
ALLSEL,ALL
F,head_pilot,FX,fx_total
F,head_pilot,FY,fy_total
ALLSEL,ALL
OUTRES,ALL,LAST
SOLVE
FINISH

/POST1
SET,LAST
*GET,uy_head_mpc,NODE,head_pilot,U,Y
*GET,ux_head_mpc,NODE,head_pilot,U,X

ETABLE,mpc_fx,SMISC,1               ! ETABLE,SMISC,1: fuerza axial local N del conector MPC184.
ETABLE,mpc_my,SMISC,2               ! ETABLE,SMISC,2: momento local My.
ETABLE,mpc_mz,SMISC,3               ! ETABLE,SMISC,3: momento local Mz.
ETABLE,mpc_mx,SMISC,4               ! ETABLE,SMISC,4: momento local Mx.
ETABLE,mpc_fz,SMISC,5               ! ETABLE,SMISC,5: fuerza cortante local V2 (Fz).
ETABLE,mpc_fy,SMISC,6               ! ETABLE,SMISC,6: fuerza cortante local V1 (Fy).
*GET,mpc_fx,ELEM,conn_elem,ETAB,mpc_fx
*GET,mpc_fy,ELEM,conn_elem,ETAB,mpc_fy
*GET,mpc_fz,ELEM,conn_elem,ETAB,mpc_fz
*GET,mpc_mx,ELEM,conn_elem,ETAB,mpc_mx
*GET,mpc_my,ELEM,conn_elem,ETAB,mpc_my
*GET,mpc_mz,ELEM,conn_elem,ETAB,mpc_mz

rfx_mpc=0
rfy_mpc=0
rfz_mpc=0
CMSEL,S,bottom_nodes
node_id=0
*DO,j,1,n_bottom
  node_id=NDNEXT(node_id)
  *GET,rfx_node,NODE,node_id,RF,FX
  *GET,rfy_node,NODE,node_id,RF,FY
  *GET,rfz_node,NODE,node_id,RF,FZ
  rfx_mpc=rfx_mpc+rfx_node
  rfy_mpc=rfy_mpc+rfy_node
  rfz_mpc=rfz_mpc+rfz_node
*ENDDO
ALLSEL,ALL
FINISH

ferr_mpc=SQRT((rfx_mpc+fx_total)**2+(rfy_mpc+fy_total)**2)/SQRT(fx_total**2+fy_total**2)

PARSAV,ALL,m12_state,par

/CLEAR,NOSTART                      ! CASO 2 - BEAM188 ELASTICO
/FILNAME,m12_beam,1
/TITLE,M12 - BEAM188 elastico
/UNITS,SI
PARRES,NEW,m12_state,par
select_tol=mesh_h*1E-4

/PREP7
ET,1,SOLID185
KEYOPT,1,2,0
ET,2,BEAM188                        ! BEAM188: viga 3D elastica como conector del tornillo (alternativa a MPC184).
MP,EX,1,young
MP,PRXY,1,nu
MP,EX,2,young
MP,PRXY,2,nu
SECTYPE,1,BEAM,CSOLID               ! SECTYPE: seccion tipo viga circular solida para BEAM188.
SECDATA,bolt_d                      ! SECDATA: diametro efectivo del tornillo (bolt_d).
TYPE,1
MAT,1
BLOCK,0,plate_x,0,plate_t,0,plate_z
BLOCK,0,plate_x,plate_t+interface_gap,2*plate_t+interface_gap,0,plate_z
div_x=MAX(1,NINT(plate_x/mesh_h))
div_y=MAX(1,NINT(plate_t/mesh_h))
div_z=MAX(1,NINT(plate_z/mesh_h))
MSHAPE,0,3D
MSHKEY,1
LSEL,S,LENGTH,,plate_x
LESIZE,ALL,,,div_x,,1
LSEL,S,LENGTH,,plate_t
LESIZE,ALL,,,div_y,,1
LSEL,S,LENGTH,,plate_z
LESIZE,ALL,,,div_z,,1
ALLSEL,ALL
VMESH,ALL

SELTOL,select_tol
NSEL,S,LOC,Y,0
CM,bottom_nodes,NODE
*GET,n_bottom,NODE,0,COUNT
NSEL,S,LOC,Y,2*plate_t+interface_gap
CM,top_load_nodes,NODE
*GET,n_top_load,NODE,0,COUNT
ALLSEL,ALL
SELTOL,

SELTOL,select_tol
NSEL,S,LOC,Y,2*plate_t+interface_gap
NSEL,R,LOC,X,bolt_x-head_absorb_r,bolt_x+head_absorb_r
NSEL,R,LOC,Z,bolt_z-head_absorb_r,bolt_z+head_absorb_r
CM,head_dep,NODE
*GET,n_head_dep_beam,NODE,0,COUNT
CMSEL,S,head_dep
NSEL,R,LOC,X,bolt_x-select_tol,bolt_x+select_tol
NSEL,R,LOC,Z,bolt_z-select_tol,bolt_z+select_tol
*GET,head_pilot,NODE,0,NUM,MIN
CMSEL,S,head_dep
CERIG,head_pilot,ALL                ! CERIG: acopla la cabecera al piloto antes del conector BEAM188.
ALLSEL,ALL

NSEL,S,LOC,Y,plate_t
NSEL,R,LOC,X,bolt_x-thread_absorb_r,bolt_x+thread_absorb_r
NSEL,R,LOC,Z,bolt_z-thread_absorb_r,bolt_z+thread_absorb_r
CM,thread_dep,NODE
*GET,n_thread_dep_beam,NODE,0,COUNT
CMSEL,S,thread_dep
NSEL,R,LOC,X,bolt_x-select_tol,bolt_x+select_tol
NSEL,R,LOC,Z,bolt_z-select_tol,bolt_z+select_tol
*GET,thread_pilot,NODE,0,NUM,MIN
CMSEL,S,thread_dep
CERIG,thread_pilot,ALL              ! CERIG: acopla la rosca al piloto inferior.
ALLSEL,ALL
SELTOL,

TYPE,2
REAL,1
MAT,2
SECNUM,1
E,head_pilot,thread_pilot
*GET,conn_elem_beam,ELEM,0,NUM,MAX
ALLSEL,ALL
FINISH

/SOLU
ANTYPE,STATIC
NLGEOM,OFF
CMSEL,S,bottom_nodes
D,ALL,ALL,0
ALLSEL,ALL
F,head_pilot,FX,fx_total
F,head_pilot,FY,fy_total
ALLSEL,ALL
OUTRES,ALL,LAST
SOLVE
FINISH

/POST1
SET,LAST
*GET,uy_head_beam,NODE,head_pilot,U,Y
*GET,ux_head_beam,NODE,head_pilot,U,X

ETABLE,beam_fx,SMISC,1              ! ETABLE,SMISC,1: fuerza axial N del conector BEAM188.
ETABLE,beam_my,SMISC,2              ! ETABLE,SMISC,2: momento My en el conector.
ETABLE,beam_mz,SMISC,3              ! ETABLE,SMISC,3: momento Mz en el conector.
ETABLE,beam_mx,SMISC,4              ! ETABLE,SMISC,4: momento Mx en el conector.
ETABLE,beam_fz,SMISC,5              ! ETABLE,SMISC,5: cortante V2 (Fz).
ETABLE,beam_fy,SMISC,6              ! ETABLE,SMISC,6: cortante V1 (Fy).
*GET,beam_fx,ELEM,conn_elem_beam,ETAB,beam_fx
*GET,beam_fy,ELEM,conn_elem_beam,ETAB,beam_fy
*GET,beam_fz,ELEM,conn_elem_beam,ETAB,beam_fz
*GET,beam_mx,ELEM,conn_elem_beam,ETAB,beam_mx
*GET,beam_my,ELEM,conn_elem_beam,ETAB,beam_my
*GET,beam_mz,ELEM,conn_elem_beam,ETAB,beam_mz

rfx_beam=0
rfy_beam=0
rfz_beam=0
CMSEL,S,bottom_nodes
node_id=0
*DO,j,1,n_bottom
  node_id=NDNEXT(node_id)
  *GET,rfx_node,NODE,node_id,RF,FX
  *GET,rfy_node,NODE,node_id,RF,FY
  *GET,rfz_node,NODE,node_id,RF,FZ
  rfx_beam=rfx_beam+rfx_node
  rfy_beam=rfy_beam+rfy_node
  rfz_beam=rfz_beam+rfz_node
*ENDDO
ALLSEL,ALL
FINISH

ferr_beam=SQRT((rfx_beam+fx_total)**2+(rfy_beam+fy_total)**2)/SQRT(fx_total**2+fy_total**2)
fx_diff=ABS(mpc_fx-beam_fx)/MAX(ABS(mpc_fx),ABS(beam_fx),1)
fy_diff=ABS(mpc_fy-beam_fy)/MAX(ABS(mpc_fy),ABS(beam_fy),1)
axial_ref=ABS(fy_total)
shear_ref=ABS(fx_total)
axial_mpc_err=ABS(ABS(mpc_fx)-axial_ref)/axial_ref
shear_mpc_err=ABS(ABS(mpc_fy)-shear_ref)/shear_ref
axial_beam_err=ABS(ABS(beam_fx)-axial_ref)/axial_ref
shear_beam_err=ABS(ABS(beam_fy)-shear_ref)/shear_ref
uy_diff=ABS(uy_head_mpc-uy_head_beam)/MAX(ABS(uy_head_mpc),ABS(uy_head_beam),1E-12)

passes=1                            ! AUDITORIA Y CONTRATO
*IF,n_nodes,LT,200,THEN
  passes=0
*ENDIF
*IF,n_solid,LT,80,THEN
  passes=0
*ENDIF
*IF,n_head_dep,LT,min_head_dep,THEN
  passes=0
*ENDIF
*IF,n_thread_dep,LT,min_thread_dep,THEN
  passes=0
*ENDIF
*IF,n_head_dep_beam,NE,n_head_dep,THEN
  passes=0
*ENDIF
*IF,n_thread_dep_beam,NE,n_thread_dep,THEN
  passes=0
*ENDIF
*IF,ferr_mpc,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,ferr_beam,GE,equilibrium_tol,THEN
  passes=0
*ENDIF
*IF,axial_mpc_err,GT,force_match_tol,THEN
  passes=0
*ENDIF
*IF,shear_mpc_err,GT,force_match_tol,THEN
  passes=0
*ENDIF
*IF,axial_beam_err,GT,force_match_tol,THEN
  passes=0
*ENDIF
*IF,shear_beam_err,GT,force_match_tol,THEN
  passes=0
*ENDIF
*IF,fx_diff,GT,force_match_tol,THEN
  passes=0
*ENDIF
*IF,fy_diff,GT,force_match_tol,THEN
  passes=0
*ENDIF

*CFOPEN,m12_connector_audit,csv
*VWRITE
('case,connector,n_nodes,n_solid,n_head_dep,n_thread_dep,N_local_N,V1_local_N,V2_local_N,mx_Nm,my_Nm,mz_Nm,uy_head_m,ux_head_m,force_error,passes')
*VWRITE,'mpc',n_nodes,n_solid,n_head_dep,n_thread_dep,mpc_fx,mpc_fy,mpc_fz,mpc_mx,mpc_my,mpc_mz,uy_head_mpc,ux_head_mpc,ferr_mpc,passes
(A4,',',F10.0,',',F10.0,',',F10.0,',',F10.0,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',F2.0)
*VWRITE,'beam',n_nodes,n_solid,n_head_dep_beam,n_thread_dep_beam,beam_fx,beam_fy,beam_fz,beam_mx,beam_my,beam_mz,uy_head_beam,ux_head_beam,ferr_beam,passes
(A4,',',F10.0,',',F10.0,',',F10.0,',',F10.0,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',F2.0)
*CFCLOS

*CFOPEN,m12_summary,csv
*VWRITE
('n_nodes,n_solid,n_head_dep,n_thread_dep,axial_mpc_err,shear_mpc_err,axial_beam_err,shear_beam_err,uy_diff_ratio,study_passes')
*VWRITE,n_nodes,n_solid,n_head_dep,n_thread_dep,axial_mpc_err,shear_mpc_err,axial_beam_err,shear_beam_err,uy_diff,passes
(F8.0,',',F8.0,',',F8.0,',',F8.0,',',F10.6,',',F10.6,',',F10.6,',',F10.6,',',F10.6,',',I1)
*CFCLOS

*CFOPEN,m12_bolt_forces,csv
*VWRITE
('connector,N_local_N,V1_local_N,V2_local_N,mx_Nm,my_Nm,mz_Nm,note')
*VWRITE,'mpc',mpc_fx,mpc_fy,mpc_fz,mpc_mx,mpc_my,mpc_mz
(A4,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',','export_for_external_code_check')
*VWRITE,'beam',beam_fx,beam_fy,beam_fz,beam_mx,beam_my,beam_mz
(A4,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',',E16.8,',','export_for_external_code_check')
*CFCLOS

*IF,passes,EQ,1,THEN
  /COM,M12 BOLT CONNECTOR STUDY PASSED
*ELSE
  /COM,M12 BOLT CONNECTOR STUDY FAILED - inspect m12_connector_audit.csv and .out
  /STATUS,PARM
*ENDIF

FINISH