"""Discrete fluid connectivity with artificial axial cuts, old and revised CAD.

Graph connectivity tests enclosure independently of a velocity/force result.
The positive-volume raw tetrahedral mesh is used only as a topology carrier.
It is not promoted to a qualified finite-volume flow mesh by this calculation.
"""
from pathlib import Path
import os,json,hashlib,sys
os.environ['OPENBLAS_NUM_THREADS']='1'
import numpy as np
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import breadth_first_order
R=Path(__file__).resolve().parents[1]
def run(path):
    raw=path.read_bytes();s=raw.decode().splitlines();ni=s.index('$Nodes');n=int(s[ni+1]);a=np.array([list(map(float,l.split())) for l in s[ni+2:ni+2+n]]);ids=a[:,0].astype(int);xyz=a[:,1:];lookup={v:i for i,v in enumerate(ids)}
    ei=s.index('$Elements');ne=int(s[ei+1]);tet=[]
    for l in s[ei+2:ei+2+ne]:
        v=list(map(int,l.split()))
        if v[1]==4:tet.append([lookup[z] for z in v[3+v[2]:]])
    tet=np.array(tet,np.int32);c=xyz[tet].mean(axis=1);n=len(tet)
    fs=np.sort(tet[:,[[0,1,2],[0,1,3],[0,2,3],[1,2,3]]].reshape(-1,3),axis=1)
    _,i,inv,cnt=np.unique(fs,axis=0,return_index=True,return_inverse=True,return_counts=True);assert max(cnt)<=2
    owner=np.repeat(np.arange(n),4);low=np.full(len(cnt),n);high=np.full(len(cnt),-1);np.minimum.at(low,inv,owner);np.maximum.at(high,inv,owner);inside=cnt==2;u=low[inside];v=high[inside]
    graph=coo_matrix((np.ones(2*len(u),np.int8),(np.r_[u,v],np.r_[v,u])),shape=(n,n)).tocsr()
    seed=np.array([.24,.076,.039]);candidate=np.flatnonzero(np.all(xyz[tet].min(axis=1)<=seed,axis=1)&np.all(xyz[tet].max(axis=1)>=seed,axis=1));root=None
    for k in candidate:
        p=xyz[tet[k]];b=np.linalg.solve((p[1:]-p[0]).T,seed-p[0])
        if min(b)>=-1e-10 and sum(b)<=1+1e-10:root=k;break
    assert root is not None,'Duct seed is outside fluid'
    cuts=[]
    for lo,hi in [(lo,.26) for lo in [.02,.06,.10,.14,.18,.22]]+[(.22,hi) for hi in [.30,.34,.38,.42,.46]]+[(.02,.46)]:
        mask=(c[:,0]>lo)&(c[:,0]<hi);index=np.flatnonzero(mask);g=graph[mask][:,mask];start=np.searchsorted(index,root);assert index[start]==root
        component=index[breadth_first_order(g,start,directed=False,return_predecessors=False)];p=c[component]
        bounded=bool(max(p[:,1])<.13 and max(p[:,2])<.13)
        reaches=bool(min(p[:,0])<lo+.006 and max(p[:,0])>hi-.006)
        cuts.append({'x_min_m':lo,'x_max_m':hi,'connected_cells':len(component),'min_centroid_x_m':float(min(p[:,0])),'max_centroid_x_m':float(max(p[:,0])),'max_centroid_y_m':float(max(p[:,1])),'max_centroid_z_m':float(max(p[:,2])),'component_remains_inside_duct_envelope':bounded,'reaches_both_axial_cuts':reaches})
    return {'mesh':str(path.relative_to(R.parent)),'mesh_sha256':hashlib.sha256(raw).hexdigest(),'seed_model_m':seed.tolist(),'cuts':cuts,'all_sampled_spans_enclosed_and_connected':all(z['component_remains_inside_duct_envelope'] and z['reaches_both_axial_cuts'] for z in cuts)}
def main():
    target=sys.argv[1] if len(sys.argv)>1 else 'p38-open-D0-U10'
    output=sys.argv[2] if len(sys.argv)>2 else 'duct-connectivity.json'
    old=run(R.parent/'cryopulse-p36/cases/p36-open-R0-U10/body-fitted.msh');new=run(R/'cases'/target/'body-fitted.msh')
    out={'physical_measurements':0,'scale':.1,'earlier_open_geometry':old,'rerouted_enclosed_candidate':new,'interpretation':'Axially clipped fluid graph components detect communication to the external domain. This is a mesh-topology and geometric-enclosure check, not a pressure seal measurement, structural proof, flow-loss estimate or driving-performance credit.'}
    assert not old['all_sampled_spans_enclosed_and_connected'],'Earlier leaking control should expose arch communication'
    assert new['all_sampled_spans_enclosed_and_connected'],'Revised geometry must retain an enclosed passage in every declared span'
    (R/'output'/output).write_text(json.dumps(out,indent=2)+'\n');print(json.dumps({'earlier_control':old['all_sampled_spans_enclosed_and_connected'],'revised_candidate':new['all_sampled_spans_enclosed_and_connected'],'sampled_spans':len(new['cuts'])}))
if __name__=='__main__':main()