"""Exploratory 2D finite-volume conduction study of E19 window installation.

Not an ISO 10211 validated solver or a product/certification calculation.
SI units. Nonoverlapping simplified sections derived from window_junction_h6.py.
Frame conductivity is calibrated against a standalone frame + calibration panel;
installation psi subtracts that SAME numerical reference, not a guessed U-window.
"""
from pathlib import Path
import argparse
import hashlib
import json
import math
import shutil
import numpy as np
from scipy import ndimage, sparse
from scipy.sparse.linalg import spsolve
from scipy.optimize import brentq

ROOT = Path(__file__).resolve().parents[3]
OUT = ROOT / 'data/site/r11/window-thermal-e19'
NAMES = ['air', 'earth', 'insulation', 'render', 'joint', 'frame', 'panel', 'support', 'wood', 'metal']
SLOPE = math.tan(math.radians(5))


def grid(a, b, step, extra=()):
    return np.unique(np.round(np.r_[np.linspace(a, b, math.ceil((b-a)/step)+1), [v for v in extra if a <= v <= b]], 9))


def geometry(kind='jamb', ret=.017, h=.002, length=1.52, support=.091,
             frame_k=.15, panel_k=.035, earth=1.2, ins=.04, strip=False):
    """x=wall depth, interior at 0; y=distance above/beside frame outer edge.

    Jamb and head are rotated to the same orientation for comparison. The sill
    keeps the actual vertical orientation. End cuts are adiabatic. Deep wall
    cells are coarser; all principal layer boundaries are inserted exactly.
    """
    xs = grid(-.03, .554, h, [0,.198,.2,.208,.21,.238,.244,.256,.262,.289,.29,.292,.295,.3,.32,.5,.503,.548,.55,.552])
    # Refine the small exposed ledge beside the inner frame, where surface
    # temperatures converge more slowly than total heat flow.
    xs = np.unique(np.r_[xs,grid(.196,.214,h/4)])
    extras = [-.08,-.06,-.028,-.025,-.023,-.02,-.017,-.005,-.003,-.002,0,.017,.04,.042,-ret-.003,-.003]
    if kind == 'reference':
        ys = grid(0,.4,h,[.04])
    elif kind == 'wall':
        ys = grid(-.2,0,h)
    else:
        ys = np.unique(np.r_[grid(-length,-.15,.02),grid(-.15,.4,h,extras),grid(-.006,.006,h/4)])
    x,y = np.meshgrid((xs[:-1]+xs[1:])/2,(ys[:-1]+ys[1:])/2)
    m = np.zeros(x.shape,dtype=np.int8)
    def put(mask, material): m[mask] = NAMES.index(material)
    def box(x0,x1,y0,y1,material): put((x>=x0)&(x<x1)&(y>=y0)&(y<y1),material)
    if kind not in ['reference']:
        top = 0 if kind=='wall' else (-ret-.003 if kind!='sill' else -.02)
        for x0,x1,mat in [(0,.2,'earth'),(.2,.3,'insulation'),(.3,.5,'earth')]:
            stop = top if kind!='sill' else {0:-.02,.2:-.06,.3:-.08}[x0]
            box(x0,x1,-length,stop,mat)
    if kind in ['jamb','head']:
        for x0,x1 in [(0,.21),(.29,.5)]:
            box(x0,x1,-ret-.003,-.003,'insulation')
            box(x0,x1,-.003,0,'render')
        box(.2,.3,-ret-.003,0,'joint')
        box(.5,.503,-ret-.003,0,'render')
    if kind=='sill':
        # Continuous load path: resolves millimetric schematic gaps/overlaps in
        # viewer. Includes insulating support beneath frame; tapes omitted.
        box(.2,.3,-.06,0,'support')
        box(-.025,.208,-.028,-.005,'wood')
        box(.5,.503,-.08,-.036,'render')
        wedge = -.017-SLOPE*(x-.29)
        put((x>=.3)&(x<.5)&(y>=-.08)&(y<wedge),'insulation')
        sheet = .017-SLOPE*(x-.29)
        if strip:
            put((x>=.32)&(x<.5)&(y>=sheet-.028)&(y<sheet-.002),'support')
        put((x>=.29)&(x<.55)&(y>=sheet-.002)&(y<sheet),'metal')
        box(.289,.292,.017,.042,'metal')
        front=.017-SLOPE*(.55-.29)
        box(.548,.552,front-.025,front,'metal')
    if kind!='wall':
        # Actual frame silhouette: 40 mm sightline, 80 mm depth. Equivalent
        # homogeneous material; NOT aluminium blocks from the 3D renderer.
        box(.21,.29,0,.04,'frame')
        box(.238,.262,.04,.4,'panel')
        if kind=='sill': box(.289,.292,.017,.042,'metal')
    if kind=='head':
        sheet=-.025+(.023/.258)*(x-.292)
        put((x>=.292)&(x<.55)&(y>=sheet)&(y<sheet+.002),'metal')
        box(.548,.552,-.002,.02,'metal')
    ks=np.array([0,earth,ins,.7,.04,frame_k,panel_k,support,.13,160.])
    return dict(xs=xs,ys=ys,m=m,k=ks[m],names=NAMES,kind=kind)


def solve(g,rsi=.13,rse=.04,ti=20.,te=0.,keep=False,ambient=None):
    """Cell-centred FVM, harmonic face conductances and Robin exposed faces.

    Air is not solved as solid. Connected air regions select interior/exterior
    boundary coefficients; no arbitrary x-position test at stepped reveals.
    """
    xs,ys,m,k=(g[n] for n in ['xs','ys','m','k'])
    solid=m>0; ids=np.full(m.shape,-1,int); n=int(solid.sum()); ids[solid]=np.arange(n)
    dx=np.broadcast_to(np.diff(xs)[None,:],m.shape)
    dy=np.broadcast_to(np.diff(ys)[:,None],m.shape)
    labels,count=ndimage.label(~solid)
    indoor=set(labels[:,0])-set([0]); outdoor=set(labels[:,-1])-set([0])
    assert not indoor & outdoor, 'Interior and exterior air communicate through geometry'
    assert set(np.unique(labels))-set([0]) == indoor|outdoor, 'Unmodelled enclosed air cavity'
    airti=np.isin(labels,list(indoor))
    diag=np.zeros(n);rhs=np.zeros(n);rows=[];cols=[];values=[];bound=[]
    for axis in [0,1]:
        aa=(slice(None,-1),slice(None)) if axis==0 else (slice(None),slice(None,-1))
        bb=(slice(1,None),slice(None)) if axis==0 else (slice(None),slice(1,None))
        distance=dy if axis==0 else dx; area=dx if axis==0 else dy
        both=solid[aa]&solid[bb]
        ia,ib=ids[aa][both],ids[bb][both]
        conduct=area[aa][both]/(distance[aa][both]/(2*k[aa][both])+distance[bb][both]/(2*k[bb][both]))
        np.add.at(diag,ia,conduct);np.add.at(diag,ib,conduct)
        rows.extend([ia,ib]);cols.extend([ib,ia]);values.extend([-conduct,-conduct])
        for ss,ar,sign in [(aa,bb,1),(bb,aa,-1)]:
            mask=solid[ss]&~solid[ar]
            idx=ids[ss][mask]; inside=airti[ar][mask]
            r=np.where(inside,rsi,rse)
            a=area[ss][mask];half=distance[ss][mask]/(2*k[ss][mask])
            c=a/(half+r);temp=np.where(inside,ti,te)
            yy,xx=np.indices(m.shape)
            px=(xs[:-1]+xs[1:])[xx[ss][mask]]/2
            py=(ys[:-1]+ys[1:])[yy[ss][mask]]/2
            if axis==0: py+=sign*distance[ss][mask]/2
            else:px+=sign*distance[ss][mask]/2
            if ambient is not None:temp=ambient(px,py,inside)
            np.add.at(diag,idx,c);np.add.at(rhs,idx,c*temp)
            bound.append((idx,c,temp,r,a,inside,m[ss][mask],px,py))
    a=sparse.coo_matrix((np.r_[np.concatenate(values),diag],(np.r_[np.concatenate(rows),np.arange(n)],np.r_[np.concatenate(cols),np.arange(n)])),shape=(n,n)).tocsr()
    t=spsolve(a,rhs)
    qi=qe=0.; surfaces=[]
    for idx,c,ta,r,area,inside,material,px,py in bound:
        flux=c*(ta-t[idx]);ts=ta-flux*r/area
        qi+=float(flux[inside].sum());qe-=float(flux[~inside].sum())
        surfaces.extend(zip(ts[inside],material[inside],px[inside],py[inside]))
    surfaces=np.array(surfaces)
    nonglass=surfaces[surfaces[:,1]!=NAMES.index('panel')]
    minface=nonglass[np.argmin(nonglass[:,0])]
    reveal=surfaces[(surfaces[:,1]!=NAMES.index('panel'))&(surfaces[:,1]!=NAMES.index('frame'))]
    result=dict(l2d=qi/(ti-te),q_in=qi,q_out=qe,balance=abs(qi-qe)/max(abs(qi),1e-12),cells=n,
                t_min=float(minface[0]),f_rsi=float((minface[0]-te)/(ti-te)),
                cold_point=dict(material=NAMES[int(minface[1])],x=float(minface[2]),y=float(minface[3])),
                t_reveal=float(reveal[:,0].min()) if len(reveal) else None)
    if keep:
        field=np.full(m.shape,np.nan);field[solid]=t
        result.update(field=field,surfaces=surfaces)
    return result


def uwall(earth=1.2,ins=.04,rsi=.13): return 1/(rsi+.4/earth+.1/ins+.04)


def analytical_checks():
    checks=[]
    wall=solve(geometry('wall'))['l2d']/.2
    checks.append(dict(name='Parede multicamada 1D: solução analítica',value=wall,reference=uwall(),
                       error=abs(wall-uwall()),ok=abs(wall-uwall())<1e-7))
    # Independent exact harmonic 2D temperature field, with matching Robin
    # boundary temperatures and adiabatic top/bottom. Exercises lateral flux.
    # T = 20-20x/L + 0.5 cos(pi*y/H) cosh(pi*(x-L/2)/H).
    L=.5;H=.4;lam=.7;k=math.pi/H
    def exact(x,y):return 20-20*x/L+.5*np.cos(k*y)*np.cosh(k*(x-L/2))
    def ambient(x,y,inside):
        derivative=-20/L+.5*k*np.cos(k*y)*np.sinh(k*(x-L/2))
        return exact(x,y)+np.where(inside,-.13,.04)*lam*derivative
    xs=np.unique(np.r_[[-.01],np.linspace(0,L,251),[L+.01]])
    ys=np.linspace(0,H,201);x,y=np.meshgrid((xs[:-1]+xs[1:])/2,(ys[:-1]+ys[1:])/2)
    m=np.where((x>0)&(x<L),1,0).astype(np.int8)
    g=dict(xs=xs,ys=ys,m=m,k=np.where(m,lam,0))
    s=solve(g,ambient=ambient,keep=True)
    err=float(np.nanmax(abs(s['field']-exact(x,y))))
    checks.append(dict(name='Campo harmónico 2D: solução analítica com fluxo lateral',max_temperature_error_K=err,ok=err<.002))
    return checks


def calibrate(uf=1.6,h=.002):
    upanel=1/(.13+.024/.035+.04)
    def residual(k):
        ref=solve(geometry('reference',h=h,frame_k=k))
        return (ref['l2d']-upanel*.36)/.04-uf
    k=brentq(residual,.005,1.,xtol=1e-9)
    reference=solve(geometry('reference',h=h,frame_k=k))
    return k,reference['l2d']


def run_case(kind,ret=.017,h=.002,uf=1.6,earth=1.2,ins=.04,support=.091,length=1.52,strip=False):
    key=(uf,h)
    if key not in CACHE:CACHE[key]=calibrate(uf,h)
    k,ref=CACHE[key]
    pars=dict(kind=kind,ret=ret,h=h,earth=earth,ins=ins,support=support,frame_k=k,length=length,strip=strip)
    heat=solve(geometry(**pars))
    # Separate real-equivalent glazing run for surface temperatures. Ug=1.1,
    # cavity/spacer not resolved; no claimed glass-edge condensation result.
    glass_k=.024/(1/1.1-.13-.04)
    surf=solve(geometry(**pars,panel_k=glass_k),rsi=.25,keep=True)
    row=dict(kind=kind,return_mm=ret*1000,mesh_mm=h*1000,uf=uf,earth_lambda=earth,ins_lambda=ins,support_lambda=support,
             frame_equivalent_lambda=k,wall_length=length,wall_u=uwall(earth,ins),strip=strip,
             reference_l2d=ref,**heat)
    row['psi']=heat['l2d']-row['wall_u']*length-ref
    row['surface']={a:b for a,b in surf.items() if a not in ['field','surfaces']}
    return row,geometry(**pars,panel_k=glass_k),surf


CACHE={}


def plot_section(name,row,g,s):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.colors import ListedColormap,BoundaryNorm
    from matplotlib.patches import Patch
    colors=['#faf9f5','#bc8c61','#e1bb54','#bbb4a4','#efe0a1','#456c6c','#a5ced3','#647d68','#997344','#8b9ea6']
    labels=['Ar','Taipa','Isolamento','Acabamento','Junta isolada','Caixilho equivalente','Vidro equivalente','Apoio isolante','Madeira','Metal']
    fig,axes=plt.subplots(1,2,figsize=(12,6.2),layout='constrained')
    fig.patch.set_facecolor('#faf9f5')
    x,y=np.meshgrid((g['xs'][:-1]+g['xs'][1:])/2,(g['ys'][:-1]+g['ys'][1:])/2)
    axes[0].pcolormesh(g['xs'],g['ys'],g['m'],cmap=ListedColormap(colors),norm=BoundaryNorm(np.arange(-.5,10.5),10),shading='flat',rasterized=True)
    heat=axes[1].pcolormesh(g['xs'],g['ys'],s['field'],cmap='coolwarm',vmin=0,vmax=20,shading='flat',rasterized=True)
    cs=axes[1].contour(x,y,np.ma.masked_invalid(s['field']),levels=[2,5,10,15,18],colors='#414745',linewidths=.65)
    axes[1].clabel(cs,fmt='%g °C',fontsize=8)
    a=s['surfaces'];r=a[(a[:,1]!=5)&(a[:,1]!=6)];cold=r[r[:,0].argmin()]
    axes[1].plot(cold[2],cold[3],'ko',ms=4)
    axes[1].annotate(f'Remate opaco: {cold[0]:.1f} °C',xy=(cold[2],cold[3]),xytext=(.025,.12),arrowprops=dict(arrowstyle='->'),fontsize=10)
    for ax in axes:
        ax.set_xlim(-.03,.555);ax.set_ylim(-.16,.15);ax.set_aspect('equal')
        ax.set_xlabel('Profundidade da parede (m) · interior ← → exterior')
        ax.set_ylabel('Distância ao bordo exterior do caixilho (m)')
        ax.set_facecolor('#faf9f5')
    axes[0].set_title('Materiais e geometria simplificada')
    axes[1].set_title('Temperaturas · interior 20 °C / exterior 0 °C')
    used=np.unique(g['m']);axes[0].legend(handles=[Patch(color=colors[i],label=labels[i]) for i in used if i],loc='upper left',fontsize=7,ncol=2)
    fig.colorbar(heat,ax=axes[1],label='°C',shrink=.65)
    title={'jamb':'Ombreira · corte rodado','head':'Verga · corte rodado','sill':'Peitoril · corte vertical'}[row['kind']]
    fig.suptitle(title+' | Estudo preliminar E.19',fontsize=15,fontweight='bold')
    fig.savefig(OUT/(name+'.png'),dpi=160)
    plt.close(fig)


def main():
    p=argparse.ArgumentParser();p.add_argument('--quick',action='store_true');args=p.parse_args()
    OUT.mkdir(parents=True,exist_ok=True)
    inputs=[ROOT/'data/site/r11/window-remates-h6/study.json',ROOT/'data/site/r11/window-remates-h6/detail/detail.json',ROOT/'tools/bim/v3/window_junction_h6.py',ROOT/'tools/bim/v3/window_remates_h6.py']
    hashes={str(p.relative_to(ROOT)):hashlib.sha256(p.read_bytes()).hexdigest() for p in inputs}
    detail=json.loads(inputs[1].read_text(encoding='utf8'))
    openings=json.loads(inputs[0].read_text(encoding='utf8'))['openings']
    assert len(openings)==15 and detail['frame_outer_m']==[.49,1.61] and detail['wall_layers_m']==[.2,.1,.2]
    results=[]
    checks=analytical_checks()
    cache_file=OUT/'run-cache.json'
    cache=json.loads(cache_file.read_text()) if cache_file.exists() else {}
    code_hash=hashlib.sha256(Path(__file__).read_bytes()).hexdigest()
    if cache.get('code_hash')!=code_hash:cache=dict(code_hash=code_hash,rows={})
    def case(name,plots=False,**kwargs):
        if name in cache['rows'] and not plots:
            row=cache['rows'][name]
        else:
            row,g,s=run_case(**kwargs)
            row['id']=name
            if plots:plot_section(name,row,g,s)
            cache['rows'][name]=row
            cache_file.write_text(json.dumps(cache),encoding='utf8')
        results.append(row)
        print(f"{name}: psi={row['psi']:.5f}; surface={row['surface']['t_reveal']:.2f} C",flush=True)
        return row
    h=.002 if args.quick else .001
    for ret in [.017,.03,.05]:
        for kind in ['jamb','head','sill']:
            if kind=='sill' and ret!=.017:continue
            case(f'{kind}-{round(ret*1000)}',kind=kind,ret=ret,h=h,plots=ret==.017)
    for uf in [1.,2.5]:
        for kind in ['jamb','sill']:
            case(f'{kind}-uf-{uf}',kind=kind,h=h,uf=uf)
    case('sill-mineral-support',kind='sill',h=h,support=1.2,plots=True)
    case('sill-through-pad',kind='sill',h=h,strip=True)
    case('jamb-lambda-plus25',kind='jamb',h=h,earth=1.5,ins=.05)
    case('jamb-no-return',kind='jamb',h=h,ret=0)
    for kind in ['jamb','head','sill']:
        coarse=case(f'{kind}-coarse',kind=kind,h=.002)
        fine=next(r for r in results if r['id']==f'{kind}-17')
        if not args.quick:
            finer=case(f'{kind}-refined',kind=kind,h=.0005)
            checks.append(dict(name=f'{kind}: malha 1 → 0.5 mm (local 0.25 → 0.125 mm)',
                delta_psi_W_mK=abs(finer['psi']-fine['psi']),delta_surface_K=abs(finer['surface']['t_reveal']-fine['surface']['t_reveal']),
                ok=abs(finer['psi']-fine['psi'])<.0002 and abs(finer['surface']['t_reveal']-fine['surface']['t_reveal'])<.1))
    long=case('jamb-long-domain',kind='jamb',h=h,length=2.02)
    base=next(r for r in results if r['id']=='jamb-17')
    checks.append(dict(name='Domínio de parede 1.52 → 2.02 m',delta_psi=abs(long['psi']-base['psi']),ok=abs(long['psi']-base['psi'])<.0001))
    balance=max(max(r['balance'],r['surface']['balance']) for r in results)
    checks.append(dict(name='Balanço de energia de todos os cenários',max_relative_residual=balance,ok=balance<1e-7))
    checks.append(dict(name='Fontes e variantes E19 preservadas',ok=all(hashlib.sha256(p.read_bytes()).hexdigest()==hashes[str(p.relative_to(ROOT))] for p in inputs)))
    checks.append(dict(name='15 janelas iguais de 49 × 161 cm e parede 20/10/20 cm confirmadas',ok=True))
    payload=dict(date='2026-10-27',status='Preliminary local study; no geometry modification or publication',
                 base_revision='r11-e19-remates-janelas-2026-10-27',inputs=hashes,results=results,checks=checks,
                 lengths_m=dict(jamb=48.3,head=7.35,sill=7.35),fixing_allowance_W_mK=.02,
                 python_libraries=dict(numpy=np.__version__,scipy=__import__('scipy').__version__))
    (OUT/'results.json').write_text(json.dumps(payload,ensure_ascii=False,indent=2),encoding='utf8')
    shutil.copyfile(Path(__file__).with_suffix('.html'),OUT/'index.html')
    assert all(c['ok'] for c in checks),checks
    print(f'{len(results)} cases; {len(checks)} checks passed. {OUT}',flush=True)


if __name__=='__main__':main()