# Grouper, round 3: the skin's mesh (exec'd after gr_skin_relax.py; same namespace): the collar's warp, the levels (the band's
# recess, the grooves' floors, the raised window frame), the faces' materials, the vertex fields of the surfaces pass (strip_v /
# strip_w, seam_chine / seam_belly, the analytic normals an_x/y/z), the skin's UV chart (UV_Grid) and the stern's duct rings.

# ================================================================== 11. the collar's square onto the bore's circle (after the relaxation)
_cq = [v for v in set(sum(([CV[i][k] for k in range(M)] for i in (I_CA, I_CB)), [])) if True]
Y_SQC = 0.5 * (Yi(STX[I_CA]) + Yi(STX[I_CB]))
A_SQ_Y = 0.5 * abs(Yi(STX[I_CA]) - Yi(STX[I_CB]))


def warp_collar_pts():
    """the grid's square round the clamp axis (stations I_CA..I_CB, columns J_HC..M-J_HC) mapped onto the circle r 0.512 in plan;
    the grid round it follows within 1.6 times the square (under the collar's skirt, r 0.695); z back on the deck's surface"""
    n = 0
    for v in range(len(SKV)):
        p = SKV[v]
        if p.z < 0.5 or abs(p.y - Y_SQC) > 1.7 * A_SQ_Y + 0.05:
            continue
        u, w_ = p.x / A_SQ, (p.y - Y_SQC) / A_SQ_Y
        rho = max(abs(u), abs(w_))
        if rho < 1e-6 or rho > 1.6:
            continue
        dx, dy = p.x, p.y - Y_SQC
        rr = math.hypot(dx, dy)
        wt = smoothstep(1.0, 1.6, rho)
        s_ = (1 - wt) * R_BORE_HOLE * rho / rr + wt
        nx, ny = dx * s_, Y_SQC + dy * s_
        SKV[v] = Vector((nx, ny, z_top(nx, ny)))
        n += 1
    return n


REPORT["collar_square"] = {"stations": [I_CA, I_CB], "half_y": round(A_SQ_Y, 4), "centre_y": round(Y_SQC, 4),
                           "warped_vertices": warp_collar_pts()}

# ================================================================== 12. levels
def floor_level(x):
    """the upper groove's floor: 50 mm under the shell beside the band, 20 mm forward of the band's front wall (the frame side
    is raised), eased between (round 2 stepped it at the wall)"""
    return LV_FLOOR_B + (LV_FLOOR_F - LV_FLOOR_B) * smoothstep(832.0, 850.0, x)


def lower_groove_depth(x):
    """the lower groove runs out on the pod's lower outer side (x 824 -> 830)"""
    return 1.0 - smoothstep(824.0, 830.0, x)


def level_of(role):
    """(level, zone) of a skin vertex: zone shell / band / floor / frame (a face whose corners share a zone is that surface; else
    it is a step's wall)"""
    kind = role[0]
    if kind == "hole":
        return {"door": (LV_BAND, "band"), "sg": (LV_FRAME, "frame"), "fg": (LV_FRAME, "frame")}[role[1]]
    j, x = role[1], role[2]
    if x < X_W - 0.05:                                          # the band zone (aft of the band's front wall, the wall included)
        if j <= J_LR[0]:
            return 0.0, "shell"
        if j in (J_LR[1], J_LR[2]):
            d = lower_groove_depth(x)
            return LV_BAND + (LV_FLOOR_B - LV_BAND) * d, ("floor" if d > 0.5 else "band")
        if j <= J_UR[0]:
            return LV_BAND, "band"
        if j in (J_UR[1], J_UR[2]):
            return floor_level(x), "floor"
        return 0.0, "shell"
    # forward of it (the transition rings, C and the cap): the frame band between FRAME_BOT and the upper groove
    if j <= J_FB[0]:
        return 0.0, "shell"
    if j <= J_UR[0]:
        return LV_FRAME, "frame"
    if j in (J_UR[1], J_UR[2]):
        return floor_level(x), "floor"
    return 0.0, "shell"


USED = set()
for fi, f in enumerate(SKF):
    if fi not in DEAD_F:
        USED.update(f)
_LZ = [level_of(SKROLE[v]) if v in USED else (0.0, "shell") for v in range(len(SKV))]
SKLV = [lz[0] for lz in _LZ]
SKZ = [lz[1] for lz in _LZ]
SKN = [proj(SKV[v])[1] if v in USED else Vector((0, 0, 1)) for v in range(len(SKV))]
for v in range(len(SKV)):                                    # the slaves' normals: mirrors of their twins' (exact symmetry)
    tw = SKTW[v]
    if tw >= 0 and v in USED:
        n_ = SKN[tw]
        SKN[v] = Vector((-n_.x, n_.y, n_.z))


def offset_dir(v):
    """the direction a vertex's level moves it: on the planar stations the surface normal within the station's plane, so that a
    station stays one plane at every level (the stern's boat-tail faces aft: a 50 mm floor along the true normal moved its
    vertices 4 cm forward, past the next station 2.2 cm away, and folded the grooves' ends); the true normal from x 826 to the
    transition rings"""
    n_ = SKN[v]
    w_ = 1.0 - smoothstep(826.0, X_W, ximg(SKV[v].y))
    if w_ <= 0.0:
        return n_
    q = Vector((n_.x, n_.y * (1.0 - w_), n_.z))
    return q.normalized() if q.length > 1e-6 else n_


SKO = [offset_dir(v) for v in range(len(SKV))]
SKP = [SKV[v] + SKO[v] * SKLV[v] for v in range(len(SKV))]


# ================================================================== 13. fields for the surfaces pass
def seam_chine(p):
    x = ximg(p.y)
    if not (545.0 <= x <= 632.0):
        return 0.0
    zc = Zi(lin_px(CHINE_SIDE, x))
    d = p.z - zc
    return 5.0 + d if abs(d) < 0.3 and abs(p.x) > 0.2 else 0.0


def seam_belly(p):
    d = p.y - Yi(BELLY_SEAM_X)
    return 5.0 + d if abs(d) < 0.3 and p.z < -0.9 else 0.0


def lin_px(tab, x, k=1):
    if x <= tab[0][0]:
        return tab[0][k]
    if x >= tab[-1][0]:
        return tab[-1][k]
    for a, b in zip(tab, tab[1:]):
        if a[0] <= x <= b[0]:
            t = (x - a[0]) / (b[0] - a[0])
            return a[k] + (b[k] - a[k]) * t
    return tab[-1][k]


# the grooves' floors: strip_v 0 / 1 on the floor's two lines, strip_w the floor's width there (its partner across the floor:
# the next vertex along the station (the planar and transition rings) or along the cap's v-line)
STRIPV = {}


def _floor_pairs():
    pairs = []
    for i in range(I_C + 1):
        for (j1, j2, lower) in ((J_UR[1], J_UR[2], False), (J_LR[1], J_LR[2], True)):
            for k1, k2 in ((j1, j2), (M - j1, M - j2)):
                a_, b_ = CV[i][k1], CV[i][k2]
                if a_ is None or b_ is None:
                    continue
                if lower and lower_groove_depth(ximg(SKV[a_].y)) < 0.5:
                    continue
                pairs.append((a_, b_))
    for a in range(NA + 1):
        a_, b_ = FV[a][J_UR[1] - J_P], FV[a][J_UR[2] - J_P]
        pairs.append((a_, b_))
    return pairs


for a_, b_ in _floor_pairs():
    w_ = (SKP[a_] - SKP[b_]).length
    STRIPV[a_] = (0.0, w_)
    STRIPV[b_] = (1.0, w_)

# ================================================================== 14. the skin's UV chart (UV_Grid, metres): the cylinder as one island (u
# round the section from the keel, port +u / starboard -u, v along the hull), the cap as a second (an azimuthal equidistant
# projection about the nose's forward axis, seen from inside the nose: injective there, no folds)
_VY = np.array(sorted(set(round(Yi(x), 6) for x in STX)))
_vcum = np.concatenate([[0.0], np.cumsum(np.abs(np.diff(_VY)))])


def _v_of(y):
    """v along the hull: the distance in y from the most forward planar station, continued (not clamped) past the planar
    stations: the transition rings and the side windows' rings lie forward of them (clamped, they all mapped onto one row of
    texels, and the bakes smeared that row over the nose's frame band: streaks and a cyan line along the ring C)"""
    return float(y - _VY[0])


def uv_cyl(p):
    S_ = sec_at(p.y)                                        # (cached per station; transition vertices have their own y)
    u = S_.nearest_arc(np.array((abs(p.x), p.z)))
    return u, _v_of(p.y)


_SEC_UV = {}


def uv_cyl_fast(p):
    key = round(p.y, 4)
    if key not in _SEC_UV:
        _SEC_UV[key] = sec_at(key)
    S_ = _SEC_UV[key]
    return S_.nearest_arc(np.array((abs(p.x), p.z))), _v_of(p.y)


R_CAP_UV = 0.70


def uv_cap(p):
    d = (p - O_NOSE).normalized()
    th = math.acos(max(-1.0, min(1.0, -d.y)))              # from the forward axis
    ph = math.atan2(d.z, d.x)
    return (R_CAP_UV * th * math.cos(ph) + 6.0, R_CAP_UV * th * math.sin(ph))   # offset: its own island beside the cylinder's


CAP_V = set()
for fi, tg in enumerate(SKFT):
    if fi in DEAD_F:
        continue
    if tg[0] == "cap" or (tg[0] == "ring" and tg[1] == "fg"):
        CAP_V.update(SKF[fi])

# ================================================================== 15. the hull mesh
hull = MB("SM_Grouper_Hull", [GUN, BAND, FRAME, STRIP, DUCTG, RECESS])
MI = {"shell": hull.mi(GUN), "band": hull.mi(BAND), "frame": hull.mi(FRAME), "floor": hull.mi(STRIP), "duct": hull.mi(DUCTG),
      "recess": hull.mi(RECESS)}
bm = bmesh.new()
lay_n = [bm.verts.layers.float.new(n) for n in ("an_x", "an_y", "an_z")]
lay_sc = bm.verts.layers.float.new("seam_chine")
lay_sb = bm.verts.layers.float.new("seam_belly")
lay_s = bm.faces.layers.int.new("on_skin")
lay_stv = bm.verts.layers.float.new("strip_v")
lay_stw = bm.verts.layers.float.new("strip_w")
lay_uvg = bm.faces.layers.int.new("uv_grid")
lay_uv = bm.loops.layers.uv.new("UV_Grid")
BV = {}


def bv(v):
    if v not in BV:
        p = SKP[v]
        b = bm.verts.new(p)
        for c, l in zip(SKN[v], lay_n):
            b[l] = c
        b[lay_sc] = seam_chine(p)
        b[lay_sb] = seam_belly(p)
        sv, sw = STRIPV.get(v, (-1.0, 0.0))
        b[lay_stv], b[lay_stw] = sv, sw
        BV[v] = b
    return BV[v]


ZRANK = {"floor": 0, "band": 1, "shell": 2, "frame": 3}


_UVC = {}
for fi, f in enumerate(SKF):
    if fi in DEAD_F:
        continue
    face = bm.faces.new([bv(v) for v in f])
    cap = all(v in CAP_V for v in f)
    sgn = 1.0 if sum(SKV[v].x for v in f) >= 0 else -1.0
    for l_, v in zip(face.loops, f):
        if cap:
            l_[lay_uv].uv = uv_cap(SKP[v])
        else:
            if v not in _UVC:
                _UVC[v] = uv_cyl_fast(SKV[v])
            u, w_ = _UVC[v]
            l_[lay_uv].uv = (sgn * u, w_)
    face[lay_uvg] = 1
    zs = set(SKZ[v] for v in f)
    if len(zs) == 1:
        role = zs.pop()
        face[lay_s] = 1
    else:                                                   # a step's wall: the higher side's material
        role = max(zs, key=lambda z: ZRANK[z])
        if role == "shell" and "band" in zs and "floor" in zs:
            role = "band"
        face[lay_s] = 0
    face.material_index = MI[role]

# ---- the stern: the duct's rings continue the grid (the strip's groove round the hull, the shroud, its rounded lip, the
# glowing inner wall, the stern plate). Angles of the last ring's columns about the duct's axis.
DUCT_C = Vector((0.0, 0.0, Zi(DUCT["y_c"])))
R_OUT = Lp(DUCT["r_out"])
Y_LIP = Yi(DUCT["x_aft"])
R_CORNER = Lp(DUCT["aft_corner_r"])
R_IN = Lp(DUCT["r_inner"]) * DUCT["r_out"] / 41.7
R_THROAT = Lp(DUCT["r_throat"]) * DUCT["r_out"] / 41.7
Y_PLATE = Yi(541.0) + 0.06
_last = [CV[0][k] for k in range(M)]
ANG = [math.atan2(SKP[v].z - DUCT_C.z, SKP[v].x) for v in _last]
R_FWD = R_OUT - 0.014
Y_STEP = Y_TAIL + 0.66 * (Y_LIP - Y_TAIL)
DUCT_RINGS = [
    (R_OUT - 0.025, Y_TAIL + 0.004, "shell"),
    (R_OUT - 0.025, Y_TAIL + 0.034, "floor"),
    (R_FWD, Y_TAIL + 0.038, "shell"),
    (R_FWD, Y_STEP - 0.010, "shell"),
    (R_FWD + 0.003, Y_STEP - 0.006, "shell"),
    (R_OUT - 0.003, Y_STEP + 0.006, "frame"),
    (R_OUT, Y_STEP + 0.010, "frame"),
    (R_OUT, Y_LIP - R_CORNER, "shell")]
for q in range(1, 5):
    a = math.pi / 2 * q / 4
    DUCT_RINGS.append((R_OUT - R_CORNER + R_CORNER * math.cos(a), Y_LIP - R_CORNER + R_CORNER * math.sin(a), "shell"))
DUCT_RINGS += [(R_IN + 0.035, Y_LIP, "frame"), (R_IN + 0.008, Y_LIP - 0.010, "frame"), (R_IN, Y_LIP - 0.035, "frame"),
               (R_THROAT, Yi(527.0), "duct"), (R_THROAT + 0.004, Yi(536.0), "duct"), (R_IN - 0.02, Y_PLATE + 0.02, "duct"),
               (R_IN - 0.045, Y_PLATE, "duct"), (Lp(DUCT["r_hub"]) + 0.03, Y_PLATE, "frame")]
prev = [bv(v) for v in _last]
R_LAST = [math.hypot(SKP[v].x, SKP[v].z - DUCT_C.z) for v in _last]
for q_, (r_, y_, role) in enumerate(DUCT_RINGS):
    if q_ == 0:
        cur = [bm.verts.new(Vector((rr * math.cos(ANG[k]), y_, DUCT_C.z + rr * math.sin(ANG[k]))))
               for k, rr in enumerate([R_LAST[k] - 0.004 for k in range(M)])]
    else:
        cur = [bm.verts.new(Vector((r_ * math.cos(ANG[k]), y_, DUCT_C.z + r_ * math.sin(ANG[k])))) for k in range(M)]
    if q_ in (0, 1):
        for v_ in cur:
            v_[lay_stv], v_[lay_stw] = float(q_), 0.030
    for k in range(M):
        f = bm.faces.new((prev[k], prev[(k + 1) % M], cur[(k + 1) % M], cur[k]))
        f.material_index = MI[role]
        f[lay_s] = 0
    prev = cur
for f in coons_cap(bm, prev, MI["frame"]):
    f[lay_s] = 0
bm.normal_update()
bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
top_f = max(bm.faces, key=lambda f: f.calc_center_median().z)
if top_f.normal.z < 0:
    bmesh.ops.reverse_faces(bm, faces=bm.faces[:])
bm.normal_update()
hull.merge(bm, None)
REPORT["hull_columns"] = M
REPORT["hull_stations"] = I_C + 1

# ---- for the later stages: the skin's points with their levels (depth under the actual skin), their normals
SKIN_PTS = [SKV[v] for v in sorted(USED)]
SKIN_LV = [SKLV[v] for v in sorted(USED)]
SKIN_NRM = [SKN[v] for v in sorted(USED)]
print("SKIN mesh", len(hull.bm.verts), len(hull.bm.faces), round(time.time() - SKIN_T0, 1), flush=True)


# ================================================================== 16. the skin's structure for the topology check (build/checks/skin_struct.json)
def _line_ids(jcol, i_from=0, i_to=None, through_cap=True):
    """a design line's vertices in order: the port column jcol over rings i_from..i_to, round the nose along the cap's u-line,
    back along the starboard column (when it runs into the cap)"""
    i_to = I_C if i_to is None else i_to
    port = [CV[i][jcol] for i in range(i_from, i_to + 1)]
    if not through_cap or i_to < I_C or not (J_P <= jcol <= J_D):
        stbd = [CV[i][M - jcol] for i in range(i_from, i_to + 1)] if 0 < jcol < J else []
        return port, stbd
    b = jcol - J_P
    cap = [FV[a][b] for a in range(NA - 1, 0, -1)]
    stbd = [CV[i][M - jcol] for i in range(I_C, i_from - 1, -1)]
    return port + cap + stbd, []


def _last_station_le(x):
    return max(i for i, xx in enumerate(STX) if xx <= x + 1e-6)


SKIN_STRUCT = {"lines": {}, "rings": [], "faces": [], "ftag": [], "pos": [], "nrm": [], "zone": [], "twin": [], "role": []}
for k in range(4):
    SKIN_STRUCT["lines"]["UR%d" % k], _ = _line_ids(J_UR[k])
for k in range(4):
    a_, b_ = _line_ids(J_LR[k], 0, _last_station_le(X_LR_END), False)
    SKIN_STRUCT["lines"]["LR%d_L" % k], SKIN_STRUCT["lines"]["LR%d_R" % k] = a_, b_
for k in range(2):
    SKIN_STRUCT["lines"]["FB%d" % k], _ = _line_ids(J_FB[k], I_BF)
SKIN_STRUCT["lines"]["stern_ring"] = [CV[0][k] for k in range(M)] + [CV[0][0]]
SKIN_STRUCT["rings"] = [[CV[i][k] for k in range(M)] for i in range(I_C + 1)]
SKIN_STRUCT["cap"] = [[FV[a][b] for b in range(N_SIDE + 1)] for a in range(NA + 1)]
SKIN_STRUCT["planar_rings"] = I_W + 1
for fi, f in enumerate(SKF):
    if fi in DEAD_F:
        continue
    SKIN_STRUCT["faces"].append(list(f))
    zs = set(SKZ[v] for v in f)
    SKIN_STRUCT["ftag"].append([SKFT[fi][0], "wall" if len(zs) > 1 else "flat"])
SKIN_STRUCT["pos"] = [[round(c, 6) for c in SKP[v]] for v in range(len(SKP))]
SKIN_STRUCT["base"] = [[round(c, 6) for c in SKV[v]] for v in range(len(SKV))]
SKIN_STRUCT["nrm"] = [[round(c, 5) for c in SKN[v]] for v in range(len(SKN))]
SKIN_STRUCT["zone"] = SKZ
SKIN_STRUCT["twin"] = SKTW
SKIN_STRUCT["kind"] = [c[0] for c in SKC]
SKIN_STRUCT["meta"] = {"M": M, "J": J, "J_UR": list(J_UR), "J_LR": list(J_LR), "J_FB": list(J_FB), "J_D": J_D, "J_P": J_P,
                       "I_W": I_W, "I_C": I_C, "I_BF": I_BF, "STX": STX}
os.makedirs(os.path.join(HERE, "checks"), exist_ok=True)