"""4.2-basalt: generate real structured 3D C3D8 thermal meshes for CalculiX.
Dimensions m, temperatures K, conductivity W/m/K. All axial faces except end
planes and all cylindrical sides insulated. This is a component boundary case.
"""
from pathlib import Path
import math,json,platform,subprocess,re,time,sys
P=Path(__file__).resolve().parent;D=json.loads((P/'design.json').read_text());O=P/'results'/'neck-fea';O.mkdir(exist_ok=True,parents=True)
r=D['neck_bore_radius_mm']/1000;L=D['neck_length_mm']/1000

def case(label,nz,nt,nr,ccx):
 folder=O/label;folder.mkdir(exist_ok=True);text=['*HEADING',f'MHeart {D["version"]} {label} component thermal case','*NODE,NSET=ALLN']
 def idx(i,j,k):return (i*(nr+1)+j)*nt+(k%nt)+1
 for i in range(nz+1):
  z=i*L/nz;t=(D['neck_wall_hot_mm']+(D['neck_wall_cold_mm']-D['neck_wall_hot_mm'])*i/nz)/1000
  for j in range(nr+1):
   rad=r+j*t/nr
   for k in range(nt):
    a=2*math.pi*k/nt;text.append(f'{idx(i,j,k)},{rad*math.cos(a):.10g},{rad*math.sin(a):.10g},{z:.10g}')
 text+=['*ELEMENT,TYPE=C3D8,ELSET=WALL'];n=0
 for i in range(nz):
  for j in range(nr):
   for k in range(nt):
    n+=1;nodes=[idx(i,j,k),idx(i,j+1,k),idx(i,j+1,k+1),idx(i,j,k+1),idx(i+1,j,k),idx(i+1,j+1,k),idx(i+1,j+1,k+1),idx(i+1,j,k+1)];text.append(','.join(map(str,[n]+nodes)))
 for name,i in [('HOT',0),('COLD',nz)]:
  text+=['*NSET,NSET='+name]
  nodes=[idx(i,j,k) for j in range(nr+1) for k in range(nt)]
  for start in range(0,len(nodes),12):text.append(','.join(map(str,nodes[start:start+12])))
 text+=['*MATERIAL,NAME=HAYNES230','*CONDUCTIVITY']
 for t,k in [(25,8.9),(100,10.4),(200,12.4),(300,14.4),(400,16.4),(500,18.4),(600,20.4),(700,22.4),(800,24.4)]:text.append(f'{k},{t+273.15}')
 text+=['*SOLID SECTION,ELSET=WALL,MATERIAL=HAYNES230','*INITIAL CONDITIONS,TYPE=TEMPERATURE','ALLN,313.15','*STEP,INC=100','*HEAT TRANSFER,STEADY STATE','1,1','*BOUNDARY','HOT,11,11,1023.15','COLD,11,11,313.15','*NODE FILE','NT,RFL','*EL FILE','HFL','*NODE PRINT,NSET=HOT,TOTALS=YES','RFL','*NODE PRINT,NSET=COLD,TOTALS=YES','RFL','*END STEP']
 (folder/'neck.inp').write_text('\n'.join(text)+'\n');start=time.monotonic()
 with (folder/'solver.log').open('w') as log:run=subprocess.run([ccx,'-i','neck'],cwd=folder,stdout=log,stderr=subprocess.STDOUT)
 result={'label':label,'host':platform.node(),'nodes':(nz+1)*(nr+1)*nt,'elements':n,'exit_code':run.returncode,'elapsed_s':time.monotonic()-start}
 if run.returncode==0:
  dat=(folder/'neck.dat').read_text();tot=[float(v) for v in re.findall(r'total heat generation for set \w+ and time[^\n]+\s+([+-]?[\d.]+E[+-]\d+)',dat)]
  result['end_heat_W']=tot[-2:]
  if len(tot)>=2:result['relative_energy_residual']=abs(tot[-1]+tot[-2])/max(abs(tot[-2]),1)
 print(result,flush=True);return result
if __name__=='__main__':
 ccx=sys.argv[1] if len(sys.argv)>1 else 'ccx'
 rows=[case('coarse',40,32,3,ccx),case('fine',80,64,5,ccx)]
 result={'version':D['version'],'scope':'Steady component thermal FEA; NOT structural/creep/coupled combustion or engine acceptance','cases':rows}
 if all(len(r.get('end_heat_W',[]))==2 for r in rows):result['mesh_relative_difference']=abs(rows[0]['end_heat_W'][0]/rows[1]['end_heat_W'][0]-1)
 (O/'summary.json').write_text(json.dumps(result,indent=2))