# Grouper, round 3: the openings' hole rings and the relaxation of the skin (exec'd after gr_skin_mesh.py; same namespace).

# ================================================================== 8. the openings' outlines (G2: no curvature jumps)
HMIN = 0.03


def side_poly(img):
    return [(Yi(x), Zi(v)) for x, v in img]


def opening_range(poly, hmin=HMIN):
    ys = [p[0] for p in poly]
    a, b = min(ys), max(ys)
    ok = []
    for i in range(1601):
        y = a + (b - a) * i / 1600
        s_ = span_at(poly, y)
        if s_ and s_[1] - s_[0] >= hmin:
            ok.append(y)
    return min(ok), max(ok)


def rounded_rect(x0, y0, x1, y1, r, n=6):
    """(kept for the later stages: the old circular-cornered rectangle)"""
    pts = []
    for (cx, cy, a0) in ((x1 - r, y0 + r, -90), (x1 - r, y1 - r, 0), (x0 + r, y1 - r, 90), (x0 + r, y0 + r, 180)):
        for k in range(n + 1):
            a = math.radians(a0 + 90.0 * k / n)
            pts.append((cx + r * math.cos(a), cy + r * math.sin(a)))
    return pts


DOOR_PX = door_outline_px(0.0, 12)
P_DOOR = side_poly(DOOR_PX)                                       # (Y, Z), the doorway's outline (G2 corners)
_sgp = fair_xy("side_glass")
P_SG = side_poly([tuple(p) for p in _sgp[::4]])                   # (Y, Z), the side glass (the sheet's fair stroke)
DOOR_YA, DOOR_YB = opening_range(P_DOOR)
SG_YA, SG_YB = opening_range(P_SG)

# the front glass: FRONT_WINDOW on SIDE's nose profile, half width 0.443 at its top (x 890.2) to 0.539 at its bottom (x 909.4);
# its corners quarter superellipses (exponent 3, G2) 7 cm along the face. Chart: (lat, Y).
X_FG0, X_FG1 = 890.2, 909.4
Y_FG0, Y_FG1 = Yi(X_FG0), Yi(X_FG1)
FG_R = 0.075


def _crest_s(y0, y1, n=60):
    s_ = 0.0
    for k in range(n):
        ya, yb = y0 + (y1 - y0) * k / n, y0 + (y1 - y0) * (k + 1) / n
        s_ += math.hypot(yb - ya, TOPf(yb) - TOPf(ya))
    return s_


def fg_hw(y):
    """the front glass's half width at station y (superellipse corners: G2)"""
    if y < Y_FG1 - 1e-6 or y > Y_FG0 + 1e-6:
        return None
    zc = TOPf(y)
    yimg = Z0 - zc / S
    t = max(0.0, min(1.0, (yimg - FRONT_WINDOW["y_top"]) / (FRONT_WINDOW["y_bot"] - FRONT_WINDOW["y_top"])))
    hw = Lp(FRONT_WINDOW["half_top_px"]) + (Lp(FRONT_WINDOW["half_bot_px"]) - Lp(FRONT_WINDOW["half_top_px"])) * t
    d = min(_crest_s(y, Y_FG0, 12), _crest_s(Y_FG1, y, 12))
    if d < FG_R:
        u = (FG_R - d) / FG_R                                     # 1 at the end, 0 where the straight side starts
        hw -= FG_R * (1.0 - (1.0 - u ** 3) ** (1.0 / 3.0))
    return max(0.05, hw)


def fg_chart(n=40):
    yy = [Y_FG1 + (Y_FG0 - Y_FG1) * (0.5 - 0.5 * math.cos(math.pi * q / n)) for q in range(n + 1)]
    right = [(fg_hw(y), y) for y in yy]
    return right + [(-l, y) for l, y in reversed(right)]


FG_POLY = fg_chart()


# ---- the hole edges in 3D (on the base surface), dense closed polylines, offset outward under the trims' lips
def _offset_closed_2d(poly, d):
    P = np.array(poly, float)
    n = len(P)
    c = P.mean(0)
    out = []
    for i in range(n):
        t = P[(i + 1) % n] - P[i - 1]
        nrm = np.array([t[1], -t[0]]) / max(1e-12, np.hypot(*t))
        if np.dot(nrm, P[i] - c) < 0:
            nrm = -nrm
        out.append(P[i] + nrm * d)
    return out


def _dense_closed(poly, step):
    P = np.array(poly + [poly[0]], float)
    seg = np.hypot(*np.diff(P, axis=0).T)
    s_ = np.concatenate([[0.0], np.cumsum(seg)])
    n = max(16, int(s_[-1] / step))
    t = np.linspace(0, s_[-1], n, endpoint=False)
    return np.stack([np.interp(t, s_, P[:, 0]), np.interp(t, s_, P[:, 1])], 1)


class Loop3(Curve3):
    """a closed 3D polyline (the hole edges)"""

    def __init__(self, P, name=""):
        P = np.asarray(P, float)
        super().__init__(np.vstack([P, P[:1]]), name)
        self.closed_n = len(P)

    def nearest(self, q, t_hint=None, window=0.4):
        if t_hint is None:
            return super().nearest(q)
        # search a window round t_hint, wrapping
        best = None
        for sh in (0.0, -self.L, self.L):
            th = t_hint + sh
            if -window <= th <= self.L + window:
                r = super().nearest(q, min(max(th, 0.0), self.L), window)
                d = (r[0] - Vector(tuple(q))).length
                if best is None or d < best[0]:
                    best = (d, r)
        return best[1]


def side_hole_loop(poly_yz, sg, off_m, step_m=0.008):
    pts = _dense_closed(_offset_closed_2d(poly_yz, off_m), step_m)
    out = []
    for (y, z) in pts:
        p, n = on_side(sg, float(y), float(z))
        out.append(tuple(p))
    return Loop3(np.array(out))


def top_hole_loop(poly_ly, off_m, step_m=0.008):
    pts = _dense_closed(_offset_closed_2d(poly_ly, off_m), step_m)
    out = []
    for (l, y) in pts:
        p, n = on_top(float(l), float(y))
        out.append(tuple(p))
    return Loop3(np.array(out))


HOLE_OFF = {"door": 0.020, "sg": 0.011, "fg": 0.011}            # under the trims' lips (door 34 mm, windows 22 mm)
DOOR_HOLE = {sg: side_hole_loop(P_DOOR, sg, HOLE_OFF["door"]) for sg in (1, -1)}
SG_HOLE = {sg: side_hole_loop(P_SG, sg, HOLE_OFF["sg"]) for sg in (1, -1)}
FG_HOLE = top_hole_loop(FG_POLY, HOLE_OFF["fg"])

# ================================================================== 9. hole rings
DEAD_F = set()                  # removed faces (indices into SKF)


def _rect_loop_cyl(i0, i1, k0, k1):
    """the boundary of the grid block stations i0..i1 x columns k0..k1 (k increasing upward on the port side), as a loop of
    vertex ids, starting at (i0, k0), going along k first"""
    L_ = []
    for k in range(k0, k1):
        L_.append(CV[i0][k % M])
    for i in range(i0, i1):
        L_.append(CV[i][k1 % M])
    for k in range(k1, k0, -1):
        L_.append(CV[i1][k % M])
    for i in range(i1, i0, -1):
        L_.append(CV[i][k0 % M])
    corners = [0, k1 - k0, (k1 - k0) + (i1 - i0), 2 * (k1 - k0) + (i1 - i0)]
    return L_, corners


def _kill_block_cyl(i0, i1, k0, k1):
    for fi, tg in enumerate(SKFT):
        if tg[0] == "cyl" and i0 <= tg[1] < i1:
            kk = tg[2]
            if k0 <= kk < k1 or (k0 > k1 and (kk >= k0 or kk < k1)):
                DEAD_F.add(fi)


def _kill_block_cap(a0, a1, b0, b1):
    for fi, tg in enumerate(SKFT):
        if tg[0] == "cap" and a0 <= tg[1] < a1 and b0 <= tg[2] < b1:
            DEAD_F.add(fi)


def _match_loop(R, corners, H):
    """the hole edge's parameters for each rectangle vertex: each corner to the hole-edge point nearest it, the sides in between
    in proportion to their lengths along the rectangle"""
    Rp = [SKV[v] for v in R]
    n = len(R)
    tc = [H.nearest(Rp[c])[1] for c in corners]
    # keep the loop's direction: the corner parameters must increase round the loop (else reverse the parameterisation)
    def unwrap(ts):
        out = [ts[0]]
        for t in ts[1:]:
            while t < out[-1]:
                t += H.L
            out.append(t)
        return out
    fwd = unwrap(tc)
    if fwd[-1] - fwd[0] > H.L:                                   # wrong way round: go backwards
        tcb = [-t for t in tc]
        bwd = unwrap(tcb)
        ts_corner = [-t for t in bwd]
        sign = -1.0
    else:
        ts_corner = fwd
        sign = 1.0
    tpar = [0.0] * n
    for q in range(4):
        ca, cb = corners[q], corners[(q + 1) % 4] if q < 3 else n
        ta = ts_corner[q]
        tb = ts_corner[(q + 1) % 4] if q < 3 else ts_corner[0] + sign * H.L
        seg = [0.0]
        for v in range(ca, cb):
            seg.append(seg[-1] + (Rp[(v + 1) % n] - Rp[v]).length)
        for m_, v in enumerate(range(ca, cb)):
            f = seg[m_] / seg[-1] if seg[-1] > 0 else 0.0
            tpar[v] = (ta + (tb - ta) * f) % H.L
    return tpar


def make_ring(R, corners, H, n_rings, role, tw_ring=None):
    """a ring of quads from the rectangle loop R (vertex ids) to the hole edge H (Loop3); returns the new loops (ids)"""
    tpar = _match_loop(R, corners, H)
    n = len(R)
    loops = [list(R)]
    for r in range(1, n_rings + 1):
        lp = []
        for q in range(n):
            ph = H.at(tpar[q])
            if r == n_rings:
                lp.append(_nv(ph, ("crv", H, tpar[q]), role))
            else:
                p = proj(SKV[R[q]].lerp(ph, r / n_rings))[0]
                lp.append(_nv(p, ("srf",), role))
        loops.append(lp)
    for a_, b_ in zip(loops, loops[1:]):
        for q in range(n):
            SKF.append((a_[q], a_[(q + 1) % n], b_[(q + 1) % n], b_[q]))
            SKFT.append(("ring", role[1]))
    return loops


def mirror_ring(loops_port, loops_stbd_rect):
    """the starboard ring: the same structure, every new vertex a slave of its port twin"""
    out = [loops_stbd_rect]
    for lp in loops_port[1:]:
        out.append([_nv(mirror(SKV[v]), ("slave",), SKROLE[v], v) for v in lp])
    n = len(loops_stbd_rect)
    for a_, b_ in zip(out, out[1:]):
        for q in range(n):
            SKF.append((a_[q], b_[q], b_[(q + 1) % n], a_[(q + 1) % n]))
            SKFT.append(("ring", "mirror"))
    return out


def _stbd_loop(R):
    """the starboard twins of a port rectangle loop (same order)"""
    out = []
    for v in R:
        tw = [u for u in range(len(SKTW)) if SKTW[u] == v]
        out.append(tw[0] if tw else v)
    return out


_TWIN_OF = {}
for _u, _t in enumerate(SKTW):
    if _t >= 0:
        _TWIN_OF[_t] = _u


def stbd_of(v):
    return _TWIN_OF.get(v, v)


# ---- the door (both flanks): the block stations I_DA..I_DB x columns J_DB..J_UR[0]
HOLES = {}
_kill_block_cyl(I_DA, I_DB, J_DB, J_UR[0])
_kill_block_cyl(I_DA, I_DB, M - J_UR[0], M - J_DB)
R_, cor_ = _rect_loop_cyl(I_DA, I_DB, J_DB, J_UR[0])
HOLES["door_L"] = make_ring(R_, cor_, DOOR_HOLE[1], 2, ("hole", "door", 0.0))
HOLES["door_R"] = mirror_ring(HOLES["door_L"], [stbd_of(v) for v in R_])
# ---- the side glass (both flanks): the transition block I_W..I_C x columns J_FB[1]..J_UR[0]
_kill_block_cyl(I_W, I_C, J_FB[1], J_UR[0])
_kill_block_cyl(I_W, I_C, M - J_UR[0], M - J_FB[1])
R_, cor_ = _rect_loop_cyl(I_W, I_C, J_FB[1], J_UR[0])
HOLES["sg_L"] = make_ring(R_, cor_, SG_HOLE[1], 2, ("hole", "sg", 0.0))
HOLES["sg_R"] = mirror_ring(HOLES["sg_L"], [stbd_of(v) for v in R_])
# ---- the front glass: the cap's block a_G..NA-a_G x b over the frame rows (FB1 .. UR0)
A_G = 2
B_G0, B_G1 = J_FB[1] - J_P, J_UR[0] - J_P
_kill_block_cap(A_G, NA - A_G, B_G0, B_G1)
_R = []
for a in range(A_G, NA - A_G):
    _R.append(FV[a][B_G0])
for b in range(B_G0, B_G1):
    _R.append(FV[NA - A_G][b])
for a in range(NA - A_G, A_G, -1):
    _R.append(FV[a][B_G1])
for b in range(B_G1, B_G0, -1):
    _R.append(FV[A_G][b])
_cor = [0, NA - 2 * A_G, (NA - 2 * A_G) + (B_G1 - B_G0), 2 * (NA - 2 * A_G) + (B_G1 - B_G0)]
HOLES["fg"] = make_ring(_R, _cor, FG_HOLE, 2, ("hole", "fg", 0.0))
# its starboard half follows its port half: ring vertex q's twin is (nA - q) (the loop starts at the starboard bottom corner);
# the two that lie on the centre line stay on it (the nose's centre line, the hole edge)
_nA, _nB = NA - 2 * A_G, B_G1 - B_G0
_N = len(_R)
for r_, lp in enumerate(HOLES["fg"][1:], start=1):
    for q, v in enumerate(lp):
        q2 = (_nA - q) % _N
        if q2 == q or (q == _nA // 2) or (q == _nA + _nB + _nA // 2):
            if r_ < len(HOLES["fg"]) - 1:
                SKC[v] = ("crv", NOSE_C, NOSE_C.nearest(SKV[v])[1])
            continue
        if SKV[v].x < 0:
            SKC[v] = ("slave",)
            SKTW[v] = lp[q2]
# ---- the collar: the deck block over the clamp axis maps onto the bore's circle (r 0.512) in plan; the grid round it follows
# within 1.6 times the square (hidden under the collar's skirt, r 0.695)
_kill_block_cyl(I_CA, I_CB, J_HC, M - J_HC)
REPORT["skin_holes"] = {k: len(v[-1]) for k, v in HOLES.items()}
print("SKIN holes", REPORT["skin_holes"], round(time.time() - SKIN_T0, 1), flush=True)


# ================================================================== 10. relaxation on the analytic surface
def _adjacency():
    rows, cols = [], []
    for fi, f in enumerate(SKF):
        if fi in DEAD_F:
            continue
        for q in range(4):
            a_, b_ = f[q], f[(q + 1) % 4]
            rows += [a_, b_]; cols += [b_, a_]
    E = np.unique(np.stack([rows, cols], 1), axis=0)
    return E[:, 0], E[:, 1]


def relax(n_iter=60, lam=0.5, report_every=20):
    er, ec = _adjacency()
    nv = len(SKV)
    deg = np.bincount(er, minlength=nv).astype(float)
    P = np.array([tuple(v) for v in SKV])
    kinds = [c[0] for c in SKC]
    free_idx = [v for v in range(nv) if kinds[v] in ("srf", "sec", "crv") and deg[v] > 0]
    sec_groups = {}
    for v in free_idx:
        if kinds[v] == "sec":
            sec_groups.setdefault(round(SKC[v][1], 6), []).append(v)
    crv_t = {v: SKC[v][2] for v in free_idx if kinds[v] == "crv"}
    slaves = [(v, SKTW[v]) for v in range(nv) if kinds[v] == "slave" and SKTW[v] >= 0]
    srf_idx = [v for v in free_idx if kinds[v] == "srf"]
    # (round 3) the order along each planar station: a free vertex stays between its column neighbours (a door ring or a
    # squeezed band could pull one past a pinned neighbour, a fold across the station)
    order = []
    for i_ in range(NI_P):
        row = [CV[i_][j_] for j_ in range(J + 1)]
        if any(r_ is None for r_ in row):
            continue
        order.append((round(Yi(STX[i_]), 6), row))
    for it in range(n_iter):
        acc = np.zeros_like(P)
        np.add.at(acc, er, P[ec])
        avg = acc / np.maximum(1.0, deg)[:, None]
        Q = P + lam * (avg - P)
        # free on the surface
        for v in srf_idx:
            P[v] = tuple(proj(Vector(tuple(Q[v])))[0])
        # on their station's section (the planar part): slide along the section
        for yk, vs in sec_groups.items():
            S_ = sec_at(yk)
            for v in vs:
                sg = 1.0 if Q[v, 0] >= 0 else -1.0
                a = S_.nearest_arc(np.array((abs(Q[v, 0]), Q[v, 2])))
                q = S_.at(a)
                P[v] = (sg * q[0], yk, q[1])
        # in order along the station (arc positions on the station's section, port half; 2 mm apart at least)
        for yk, row in order:
            S_ = sec_at(yk)
            arcs = [S_.nearest_arc(np.array((abs(P[v, 0]), P[v, 2]))) for v in row]
            for q in range(1, len(row) - 1):
                v = row[q]
                if kinds[v] != "sec":
                    continue
                lo_, hi_ = arcs[q - 1] + 0.002, arcs[q + 1] - 0.002
                if lo_ < hi_ and not (lo_ <= arcs[q] <= hi_):
                    a = min(max(arcs[q], lo_), hi_)
                    arcs[q] = a
                    qq = S_.at(a)
                    P[v] = (qq[0], yk, qq[1])
        # on their design line
        for v, t in crv_t.items():
            C3 = SKC[v][1]
            q, t2 = C3.nearest(Vector(tuple(Q[v])), t, 0.25)
            crv_t[v] = t2
            P[v] = tuple(q)
        for v, tw in slaves:
            P[v] = (-P[tw, 0], P[tw, 1], P[tw, 2])
        if report_every and (it % report_every == 0 or it == n_iter - 1):
            print("RELAX", it, round(float(np.abs(Q - P).max()), 5), round(time.time() - SKIN_T0, 1), flush=True)
    for v in range(nv):
        SKV[v] = Vector(tuple(P[v]))
    for v, t in crv_t.items():
        SKC[v] = (SKC[v][0], SKC[v][1], t)


relax(int(os.environ.get("GROUPER_RELAX", "60")))
REPORT["skin_relax_s"] = round(time.time() - SKIN_T0, 1)