# Grouper: helper library (exec'd first by build_grouper.py). Copied 2026-10-22 from the Mako's build/lib2.py,
# build/tender_lib.py and build/mk_topo.py (models/tender-mako, not edited there): the mesh builder MB (quads, material
# slots, merged meshes), lathe/revolve, grid_rings, foil, loft, Surface, prune_inside, Coons caps, cube spheres, ladder
# patches. Grouper additions in gr_parts.py.

# ====================================================================== from lib2.py
# Topology-first helper library (no booleans, quads, merged meshes with material slots).
import bpy, bmesh, math, os
from mathutils import Vector, Matrix, Euler


def get_coll(name, parent=None):
    c = bpy.data.collections.get(name)
    if c is None:
        c = bpy.data.collections.new(name)
        (parent or bpy.context.scene.collection).children.link(c)
    return c

def MAT(name, color=(0.5, 0.5, 0.5), metal=0.0, rough=0.5, emit=None, emit_str=0.0, alpha=1.0, coat=0.0):
    m = bpy.data.materials.get(name)
    if m is not None:
        return m
    m = bpy.data.materials.new(name)
    m.use_nodes = True
    b = m.node_tree.nodes.get("Principled BSDF")
    c = tuple(color) + (1.0,) if len(color) == 3 else tuple(color)
    b.inputs["Base Color"].default_value = c
    b.inputs["Metallic"].default_value = metal
    b.inputs["Roughness"].default_value = rough
    if "Coat Weight" in b.inputs:
        b.inputs["Coat Weight"].default_value = coat
    if emit is not None:
        b.inputs["Emission Color"].default_value = tuple(emit) + (1.0,)
        b.inputs["Emission Strength"].default_value = emit_str
    if alpha < 1.0:
        b.inputs["Alpha"].default_value = alpha
        try:
            m.surface_render_method = 'BLENDED'
        except Exception:
            pass
    m.diffuse_color = c
    return m

def clean_scene(keep=("DNR_CNC_Root", "Camera", "Light")):
    for o in list(bpy.data.objects):
        if o.name in keep or o.name.startswith(("Key_", "Fill_", "Rim_", "Top_", "Floor")):
            continue
        bpy.data.objects.remove(o)
    for me in list(bpy.data.meshes):
        if me.users == 0:
            bpy.data.meshes.remove(me)
    for cu in list(bpy.data.curves):
        if cu.users == 0:
            bpy.data.curves.remove(cu)

# ------------------------------------------------------------------ transforms
def _M(loc=(0, 0, 0), rot=(0, 0, 0), scale=(1, 1, 1)):
    S = Matrix.Diagonal(Vector(scale).to_4d())
    return Matrix.Translation(loc) @ Euler(rot).to_matrix().to_4x4() @ S

def bevel_edges(bm, edges, w, seg=1):
    if w > 0 and edges:
        bmesh.ops.bevel(bm, geom=edges, offset=w, offset_type='OFFSET', segments=seg,
                        profile=0.5, affect='EDGES')

_AX = {'+x': (0, 1), '-x': (0, -1), '+y': (1, 1), '-y': (1, -1), '+z': (2, 1), '-z': (2, -1)}
def _drop_faces(bm, drop, size):
    """delete the large flat face(s) of a box on the named sides (before bevel)."""
    fs = []
    for d in drop:
        ax, sg = _AX[d]
        for f in bm.faces:
            n = f.normal
            if abs(n[ax] - sg) < 1e-4:
                fs.append(f)
    return fs

class MB:
    """Mesh builder: accumulates primitives (quad topology) into one mesh with material slots."""
    def __init__(self, name, mats):
        self.name = name
        self.bm = bmesh.new()
        self.mats = list(mats)

    @classmethod
    def load(cls, name):
        o = bpy.data.objects[name]
        mb = cls(name, list(o.data.materials))
        mb.bm.from_mesh(o.data)
        mb._old = (o.parent, o.location.copy(), o.users_collection[0], o)
        return mb

    def dome(self, r, ctr, mat, segs=16, rings=8, rot=(0, 0, 0)):
        bm = bmesh.new()
        bmesh.ops.create_uvsphere(bm, u_segments=segs, v_segments=rings * 2, radius=r)
        low = [v for v in bm.verts if v.co.z < -1e-6]
        bmesh.ops.delete(bm, geom=low, context='VERTS')
        self.merge(bm, mat, _M(ctr, rot))

    def mi(self, m):
        if m not in self.mats:
            self.mats.append(m)
        return self.mats.index(m)

    def merge(self, bm2, mat, xf=None):
        if xf is not None:
            bm2.transform(xf)
        if mat is not None:
            mi = self.mi(mat)
            for f in bm2.faces:
                f.material_index = mi
        tmp = bpy.data.meshes.new("_t")
        bm2.to_mesh(tmp)
        bm2.free()
        self.bm.from_mesh(tmp)
        bpy.data.meshes.remove(tmp)

    # ---- primitives
    def box(self, size, ctr, mat, rot=(0, 0, 0), bevel=0.0, seg=1, drop=()):
        bm = bmesh.new()
        bmesh.ops.create_cube(bm, size=1.0)
        for v in bm.verts:
            v.co.x *= size[0]; v.co.y *= size[1]; v.co.z *= size[2]
        bm.normal_update()
        dfaces = _drop_faces(bm, drop, size)
        if bevel > 0:
            es = [e for e in bm.edges if not any(f in dfaces for f in e.link_faces)]
            bevel_edges(bm, es, bevel, seg)
            bm.normal_update()
            dfaces = _drop_faces(bm, drop, size)
        if dfaces:
            bmesh.ops.delete(bm, geom=dfaces, context='FACES_ONLY')
        self.merge(bm, mat, _M(ctr, rot))

    def cyl(self, r1, r2, h, ctr, mat, rot=(0, 0, 0), segs=24, caps='both', bevel=0.0):
        bm = bmesh.new()
        bmesh.ops.create_cone(bm, cap_ends=True, cap_tris=False, segments=segs, radius1=r1, radius2=r2, depth=h)
        bm.normal_update()
        if bevel > 0:
            es = [e for e in bm.edges if len(e.link_faces) == 2 and abs(e.verts[0].co.z - e.verts[1].co.z) < 1e-6]
            bevel_edges(bm, es, bevel, 1)
            bm.normal_update()
        rm = []
        for f in bm.faces:
            if f.normal.z > 0.999 and caps in ('bottom', 'none'):
                rm.append(f)
            if f.normal.z < -0.999 and caps in ('top', 'none'):
                rm.append(f)
        if rm:
            bmesh.ops.delete(bm, geom=rm, context='FACES_ONLY')
        self.merge(bm, mat, _M(ctr, rot))

    def sphere(self, r, ctr, mat, segs=20, rings=12, scale=(1, 1, 1), rot=(0, 0, 0)):
        bm = bmesh.new()
        bmesh.ops.create_uvsphere(bm, u_segments=segs, v_segments=rings, radius=r)
        self.merge(bm, mat, _M(ctr, rot, scale))

    def torus(self, R, r, ctr, mat, rot=(0, 0, 0), segs=40, rsegs=8):
        bm = bmesh.new()
        rows = []
        for i in range(segs):
            a = 2 * math.pi * i / segs
            rows.append([bm.verts.new(((R + r * math.cos(2 * math.pi * j / rsegs)) * math.cos(a),
                                      (R + r * math.cos(2 * math.pi * j / rsegs)) * math.sin(a),
                                      r * math.sin(2 * math.pi * j / rsegs))) for j in range(rsegs)])
        for i in range(segs):
            for j in range(rsegs):
                bm.faces.new((rows[i][j], rows[(i + 1) % segs][j], rows[(i + 1) % segs][(j + 1) % rsegs], rows[i][(j + 1) % rsegs]))
        bm.normal_update()
        self.merge(bm, mat, _M(ctr, rot))

    def quad(self, p0, p1, p2, p3, mat):
        bm = bmesh.new()
        bm.faces.new([bm.verts.new(p) for p in (p0, p1, p2, p3)])
        self.merge(bm, mat)

    def prism(self, pts, thick, ctr, mat, rot=(0, 0, 0), bevel=0.0, drop=()):
        """convex/simple polygon (u,v) in XZ plane, extruded along +Y by thick (front face at y=0)."""
        bm = bmesh.new()
        vs = [bm.verts.new((u, 0, v)) for u, v in pts]
        f = bm.faces.new(vs)
        res = bmesh.ops.extrude_face_region(bm, geom=[f])
        for v in [e for e in res['geom'] if isinstance(e, bmesh.types.BMVert)]:
            v.co.y += thick
        bm.normal_update()
        bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
        if bevel > 0:
            bevel_edges(bm, bm.edges[:], bevel, 1)
        self.merge(bm, mat, _M(ctr, rot))

    def sweep(self, pts, radii, mat, sides=8, caps=True, up=(0, 0, 1)):
        """tube along a polyline with per-point radius (quad rings, tri-fan caps)."""
        bm = bmesh.new()
        P = [Vector(p) for p in pts]
        n = len(P)
        if not isinstance(radii, (list, tuple)):
            radii = [radii] * n
        # parallel transport frame
        t0 = (P[1] - P[0]).normalized()
        ref = Vector(up)
        if abs(ref.dot(t0)) > 0.95:
            ref = Vector((1, 0, 0))
        nrm = (ref - ref.dot(t0) * t0).normalized()
        rings = []
        prev_t = t0
        for i in range(n):
            if i == 0:
                t = t0
            elif i == n - 1:
                t = (P[i] - P[i - 1]).normalized()
            else:
                t = ((P[i] - P[i - 1]).normalized() + (P[i + 1] - P[i]).normalized()).normalized()
            # rotate nrm from prev_t to t
            axis = prev_t.cross(t)
            if axis.length > 1e-7:
                ang = math.acos(max(-1, min(1, prev_t.dot(t))))
                nrm = (Matrix.Rotation(ang, 3, axis.normalized()) @ nrm)
            nrm = (nrm - nrm.dot(t) * t).normalized()
            bn = t.cross(nrm)
            prev_t = t
            ring = []
            for k in range(sides):
                a = 2 * math.pi * k / sides
                ring.append(bm.verts.new(P[i] + radii[i] * (math.cos(a) * nrm + math.sin(a) * bn)))
            rings.append(ring)
        uvl = bm.loops.layers.uv.new("UVMap")
        acc = [0.0]
        for i in range(1, n):
            acc.append(acc[-1] + (P[i] - P[i - 1]).length)
        vs = 1.0 / (2 * math.pi * (sum(radii) / n))
        for i in range(n - 1):
            for k in range(sides):
                f = bm.faces.new((rings[i][k], rings[i + 1][k], rings[i + 1][(k + 1) % sides], rings[i][(k + 1) % sides]))
                uvs = ((k / sides, acc[i] * vs), (k / sides, acc[i + 1] * vs), ((k + 1) / sides, acc[i + 1] * vs), ((k + 1) / sides, acc[i] * vs))
                for l, uv in zip(f.loops, uvs):
                    l[uvl].uv = uv
        if caps:
            for ring, c in ((rings[0], P[0]), (rings[-1], P[-1])):
                cv = bm.verts.new(c)
                for k in range(sides):
                    f = bm.faces.new((ring[k], ring[(k + 1) % sides], cv))
                    for l in f.loops:
                        l[uvl].uv = (0.5, 0.5)
        bm.normal_update()
        bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
        self.merge(bm, mat)

    def add_bm(self, bm2, mat=None, xf=None):
        self.merge(bm2, mat, xf)

    # ---- finish
    def build(self, coll, loc=(0, 0, 0), parent=None, sharp_deg=35.0, tri_ngons=True, weld=1e-5, recalc=False, origin=None):
        bm = self.bm
        bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=weld)
        bm.faces.ensure_lookup_table()
        zf = [f for f in bm.faces if f.calc_area() < 1e-10]
        self.zero_removed = len(zf)
        if zf:
            bmesh.ops.delete(bm, geom=zf, context='FACES')
        if recalc:
            bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
        # drop loose verts/edges
        loose = [v for v in bm.verts if not v.link_faces]
        if loose:
            bmesh.ops.delete(bm, geom=loose, context='VERTS')
        if tri_ngons:
            ng = [f for f in bm.faces if len(f.verts) > 4]
            if ng:
                bmesh.ops.triangulate(bm, faces=ng, quad_method='BEAUTY', ngon_method='BEAUTY')
        lim = math.radians(sharp_deg)
        bm.normal_update()
        lay_c = bm.faces.layers.int.get("crisp")          # Mako v04: faces flagged crisp (faceted parts) get all edges sharp
        for e in bm.edges:
            if len(e.link_faces) == 2:
                e.smooth = e.link_faces[0].normal.angle(e.link_faces[1].normal) < lim
                if lay_c is not None and (e.link_faces[0][lay_c] or e.link_faces[1][lay_c]):
                    e.smooth = False
            else:
                e.smooth = True
        for f in bm.faces:
            f.smooth = True
        if origin is not None:
            bmesh.ops.translate(bm, verts=bm.verts[:], vec=-Vector(origin))
            loc = origin
        me = bpy.data.meshes.new(self.name)
        bm.to_mesh(me)
        bm.free()
        for m in self.mats:
            me.materials.append(m)
        old = getattr(self, '_old', None)
        if old is not None:
            parent = parent or old[0]
            if origin is None:
                loc = old[1]
            kids = list(old[3].children)
            oldme = old[3].data
            bpy.data.objects.remove(old[3])
            bpy.data.meshes.remove(oldme)
        else:
            kids = []
        o = bpy.data.objects.new(self.name, me)
        coll.objects.link(o)
        o.location = loc
        for k in kids:
            k.parent = o
            k.matrix_parent_inverse = o.matrix_world.inverted()
        if parent is not None:
            bpy.context.view_layer.update()
            o.parent = parent
            o.matrix_parent_inverse = parent.matrix_world.inverted()
        return o

# ------------------------------------------------------------------ box/loop-cut utilities for the housing
def cube_bm(size, ctr, bevel=0.0):
    bm = bmesh.new()
    bmesh.ops.create_cube(bm, size=1.0)
    for v in bm.verts:
        v.co.x = v.co.x * size[0] + ctr[0]; v.co.y = v.co.y * size[1] + ctr[1]; v.co.z = v.co.z * size[2] + ctr[2]
    if bevel > 0:
        bevel_edges(bm, bm.edges[:], bevel, 1)
    return bm

def loopcut(bm, axis, value):
    """cut every face crossing the plane axis=value (like a loop cut through the whole mesh)."""
    n = Vector((1 if axis == 'x' else 0, 1 if axis == 'y' else 0, 1 if axis == 'z' else 0))
    bmesh.ops.bisect_plane(bm, geom=bm.verts[:] + bm.edges[:] + bm.faces[:], plane_co=n * value, plane_no=n,
                           dist=1e-6, clear_inner=False, clear_outer=False)

def faces_where(bm, nrm, plane_ax, plane_val, rng_a, rng_b, tol=1e-4):
    """faces with normal ~ nrm lying on plane (axis index, value) whose centre is inside the (a,b) ranges."""
    a_ax, b_ax = [i for i in range(3) if i != plane_ax]
    out = []
    for f in bm.faces:
        if f.normal.dot(Vector(nrm)) < 0.999:
            continue
        c = f.calc_center_median()
        if abs(c[plane_ax] - plane_val) > tol:
            continue
        if rng_a[0] < c[a_ax] < rng_a[1] and rng_b[0] < c[b_ax] < rng_b[1]:
            out.append(f)
    return out

def punch_and_extrude(bm, faces, vec):
    """delete faces and extrude the hole boundary along vec (a clean reveal tunnel)."""
    cand = set(e for f in faces for e in f.edges)
    bmesh.ops.delete(bm, geom=faces, context='FACES_ONLY')
    edges = [e for e in cand if e.is_valid and len(e.link_faces) == 1]
    res = bmesh.ops.extrude_edge_only(bm, edges=edges)
    vs = [g for g in res['geom'] if isinstance(g, bmesh.types.BMVert)]
    for v in vs:
        v.co += Vector(vec)
    bm.normal_update()
    return vs

# ------------------------------------------------------------------ audits
def audit(obj):
    me = obj.data
    bm = bmesh.new(); bm.from_mesh(me)
    tris = quads = ngons = 0
    for f in bm.faces:
        n = len(f.verts)
        if n == 3: tris += 1
        elif n == 4: quads += 1
        else: ngons += 1
    nonman = sum(1 for e in bm.edges if len(e.link_faces) != 2)
    ntri = sum(max(1, len(f.verts) - 2) for f in bm.faces)
    zero = sum(1 for f in bm.faces if f.calc_area() < 1e-9)
    loose = sum(1 for v in bm.verts if not v.link_faces)
    bm.free()
    return dict(name=obj.name, verts=len(me.vertices), tris_total=ntri, tri=tris, quad=quads, ngon=ngons,
                open_edges=nonman, zero_area=zero, loose=loose, uv=len(me.uv_layers))


def pocket(bm, face, inset, depth):
    """inset a face and sink the inner part with vertical walls (inner face is the floor)."""
    n = face.normal.copy()
    bmesh.ops.inset_individual(bm, faces=[face], thickness=inset, use_even_offset=True)
    res = bmesh.ops.extrude_face_region(bm, geom=[face])
    nf = [g for g in res['geom'] if isinstance(g, bmesh.types.BMFace)]
    nv = [g for g in res['geom'] if isinstance(g, bmesh.types.BMVert)]
    for v in nv:
        v.co -= n * depth
    bmesh.ops.delete(bm, geom=[face], context='FACES')
    bm.normal_update()
    return nf


def _mb_transform(self, M):
    bmesh.ops.transform(self.bm, matrix=M, verts=self.bm.verts[:])
MB.transform = _mb_transform

def catmull(pts, per=6):
    """Catmull-Rom resample of control points (list of 3-tuples) -> list of Vectors."""
    P = [Vector(p) for p in pts]
    P = [P[0] * 2 - P[1]] + P + [P[-1] * 2 - P[-2]]
    out = []
    for i in range(1, len(P) - 2):
        p0, p1, p2, p3 = P[i - 1], P[i], P[i + 1], P[i + 2]
        for k in range(per):
            t = k / per
            t2, t3 = t * t, t * t * t
            out.append(0.5 * ((2 * p1) + (-p0 + p2) * t + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t2 + (-p0 + 3 * p1 - 3 * p2 + p3) * t3))
    out.append(P[-2])
    return out

def sink(bm, face, depth):
    """sink a face with vertical walls (no inset); returns the new floor face."""
    n = face.normal.copy()
    res = bmesh.ops.extrude_face_region(bm, geom=[face])
    nf = [g for g in res['geom'] if isinstance(g, bmesh.types.BMFace)]
    for v in [g for g in res['geom'] if isinstance(g, bmesh.types.BMVert)]:
        v.co -= n * depth
    bmesh.ops.delete(bm, geom=[face], context='FACES')
    bm.normal_update()
    return nf[0]

# ====================================================================== from tender_lib.py
# Helpers for the Tenders (copied from tender-manta) (exec'd after lib2.py, shares its namespace).
import bpy, bmesh, math
from mathutils import Vector, Matrix
from mathutils.bvhtree import BVHTree
from mathutils import kdtree
from mathutils.interpolate import poly_3d_calc

YPOS = Matrix.Rotation(-math.pi / 2, 4, 'X')        # lathe axis (local +Z) -> world +Y
T = lambda *p: Matrix.Translation(p)


def pchip(pts):
    """monotone cubic through [(x, y)]; clamped outside the table"""
    xs = [p[0] for p in pts]; ys = [p[1] for p in pts]
    n = len(xs)
    h = [xs[i + 1] - xs[i] for i in range(n - 1)]
    d = [(ys[i + 1] - ys[i]) / h[i] for i in range(n - 1)]
    m = [d[0]] + [0.0] * (n - 2) + [d[-1]]
    for i in range(1, n - 1):
        if d[i - 1] * d[i] > 0:
            w1, w2 = 2 * h[i] + h[i - 1], h[i] + 2 * h[i - 1]
            m[i] = (w1 + w2) / (w1 / d[i - 1] + w2 / d[i])
    def f(x):
        if x <= xs[0]:
            return ys[0]
        if x >= xs[-1]:
            return ys[-1]
        lo, hi = 0, n - 1
        while hi - lo > 1:
            mid = (lo + hi) // 2
            if xs[mid] <= x:
                lo = mid
            else:
                hi = mid
        i = lo
        t = (x - xs[i]) / h[i]
        t2, t3 = t * t, t * t * t
        return ((2 * t3 - 3 * t2 + 1) * ys[i] + (t3 - 2 * t2 + t) * h[i] * m[i] + (-2 * t3 + 3 * t2) * ys[i + 1]
                + (t3 - t2) * h[i] * m[i + 1])
    return f


def cspline(pts):
    """(round 3) natural cubic spline through [(x, y)] (C2: curvature continuous, unlike pchip whose curvature jumps at
    every knot); clamped (constant) outside the table"""
    xs = [float(p[0]) for p in pts]; ys = [float(p[1]) for p in pts]
    n = len(xs)
    h = [xs[i + 1] - xs[i] for i in range(n - 1)]
    # tridiagonal system for the second derivatives (natural ends)
    a = [0.0] * n; b = [1.0] * n; c = [0.0] * n; r = [0.0] * n
    for i in range(1, n - 1):
        a[i] = h[i - 1]; b[i] = 2.0 * (h[i - 1] + h[i]); c[i] = h[i]
        r[i] = 6.0 * ((ys[i + 1] - ys[i]) / h[i] - (ys[i] - ys[i - 1]) / h[i - 1])
    for i in range(1, n):                                # Thomas algorithm
        w = a[i] / b[i - 1]
        b[i] -= w * c[i - 1]; r[i] -= w * r[i - 1]
    m = [0.0] * n
    m[-1] = r[-1] / b[-1]
    for i in range(n - 2, -1, -1):
        m[i] = (r[i] - c[i] * m[i + 1]) / b[i]

    def f(x):
        if x <= xs[0]:
            return ys[0]
        if x >= xs[-1]:
            return ys[-1]
        lo, hi = 0, n - 1
        while hi - lo > 1:
            mid = (lo + hi) // 2
            if xs[mid] <= x:
                lo = mid
            else:
                hi = mid
        i = lo
        t = x - xs[i]; hh = h[i]
        A = (m[i + 1] - m[i]) / (6.0 * hh)
        B = m[i] / 2.0
        C = (ys[i + 1] - ys[i]) / hh - hh * (2.0 * m[i] + m[i + 1]) / 6.0
        return ys[i] + t * (C + t * (B + t * A))
    return f


def smoothstep(a, b, x):
    t = max(0.0, min(1.0, (x - a) / (b - a)))
    return t * t * (3 - 2 * t)


def lathe(mb, prof, mats, segs=32, xf=None):
    """Revolve a profile [(r, a)] about local +Z (a = height). Starts and ends on the axis (closed shell, fan ends).
    Writes UV_Tile in metres (flat bands seen along the axis, walls unrolled). Same as the docking arm's."""
    bm = bmesh.new()
    uvl = bm.loops.layers.uv.new("UV_Tile")
    cols = []
    for (r, a) in prof:
        if r < 1e-9:
            cols.append([bm.verts.new((0, 0, a))])
        else:
            cols.append([bm.verts.new((r * math.cos(2 * math.pi * k / segs), r * math.sin(2 * math.pi * k / segs), a))
                         for k in range(segs)])
    ml = list(mats) if isinstance(mats, (list, tuple)) else [mats] * (len(prof) - 1)
    assert len(ml) == len(prof) - 1, (len(ml), len(prof))
    sp = [0.0]
    for i in range(len(prof) - 1):
        sp.append(sp[-1] + math.hypot(prof[i + 1][0] - prof[i][0], prof[i + 1][1] - prof[i][1]))
    for i in range(len(prof) - 1):
        A, B = cols[i], cols[i + 1]
        mi = mb.mi(ml[i])
        flat = abs(prof[i + 1][0] - prof[i][0]) >= abs(prof[i + 1][1] - prof[i][1])
        rm = (prof[i][0] + prof[i + 1][0]) / 2
        for k in range(segs):
            k2 = (k + 1) % segs
            if len(A) == 1 and len(B) == 1:
                continue
            if len(A) == 1:
                f = bm.faces.new((A[0], B[k], B[k2]))
            elif len(B) == 1:
                f = bm.faces.new((A[k], B[0], A[k2]))
            else:
                f = bm.faces.new((A[k], B[k], B[k2], A[k2]))
            f.material_index = mi
            if flat:
                for l in f.loops:
                    l[uvl].uv = (l.vert.co.x, l.vert.co.y)
            else:
                u_a, u_b = rm * 2 * math.pi * k / segs, rm * 2 * math.pi * (k + 1) / segs
                uvs = {}
                for col, s in ((A, sp[i]), (B, sp[i + 1])):
                    if len(col) == 1:
                        uvs[col[0]] = ((u_a + u_b) / 2, s)
                    else:
                        uvs[col[k]] = (u_a, s); uvs[col[k2]] = (u_b, s)
                for l in f.loops:
                    l[uvl].uv = uvs[l.vert]
    bm.normal_update()
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    mb.merge(bm, None, xf)


def revolve_loop(mb, loop, mats, segs=32, xf=None):
    """Revolve a CLOSED profile loop [(r, a)] (no point on the axis) about local +Z: a ring torus-like solid."""
    bm = bmesh.new()
    uvl = bm.loops.layers.uv.new("UV_Tile")
    cols = [[bm.verts.new((r * math.cos(2 * math.pi * k / segs), r * math.sin(2 * math.pi * k / segs), a)) for k in range(segs)]
            for (r, a) in loop]
    n = len(loop)
    ml = list(mats) if isinstance(mats, (list, tuple)) else [mats] * n
    sp = [0.0]
    for i in range(n):
        j = (i + 1) % n
        sp.append(sp[-1] + math.hypot(loop[j][0] - loop[i][0], loop[j][1] - loop[i][1]))
    for i in range(n):
        j = (i + 1) % n
        A, B = cols[i], cols[j]
        mi = mb.mi(ml[i])
        rm = (loop[i][0] + loop[j][0]) / 2
        for k in range(segs):
            k2 = (k + 1) % segs
            f = bm.faces.new((A[k], B[k], B[k2], A[k2]))
            f.material_index = mi
            u_a, u_b = rm * 2 * math.pi * k / segs, rm * 2 * math.pi * (k + 1) / segs
            uvs = {A[k]: (u_a, sp[i]), A[k2]: (u_b, sp[i]), B[k]: (u_a, sp[i + 1]), B[k2]: (u_b, sp[i + 1])}
            for l in f.loops:
                l[uvl].uv = uvs[l.vert]
    bm.normal_update()
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    mb.merge(bm, None, xf)


def grid_rings(mb, rings, mats, closed_u=True, cap0=None, cap1=None, uv=None, skip=None, expect=None, vattr=None):
    """Skin consecutive rings (lists of points, equal length) with quads. mats: one material or one per band.
    cap0/cap1: a point to fan the first/last ring to. uv: rings of (u, v) for UV_Tile (cap points take uv_cap0/1 from
    the ring's first entry). skip(i, k) -> True leaves out the face between ring i and i+1 at k (i = -1: cap0 fan,
    i = len-1: cap1 fan). expect = (i, k, direction): the faces are turned so that face (i, k) faces that way."""
    bm = bmesh.new()
    uvl = bm.loops.layers.uv.new("UV_Tile")
    lays = {name: bm.verts.layers.float.new(name) for name in (vattr or {})}   # layers first: a new layer invalidates verts
    V = [[bm.verts.new(p) for p in r] for r in rings]
    for name, vals in (vattr or {}).items():                  # per-vertex float fields (rings of values), e.g. seam distances
        lay = lays[name]
        for r, vr in zip(V, vals):
            for v, x in zip(r, vr):
                v[lay] = x
    m = len(rings[0])
    nr = len(rings)
    ml = list(mats) if isinstance(mats, (list, tuple)) else [mats] * (nr - 1)
    span = m if closed_u else m - 1
    probe = None
    for i in range(nr - 1):
        mi = mb.mi(ml[i])
        for k in range(span):
            if skip is not None and skip(i, k):
                continue
            k2 = (k + 1) % m
            ids = [(i, k), (i + 1, k), (i + 1, k2), (i, k2)]
            f = bm.faces.new([V[a][b] for a, b in ids])
            f.material_index = mi
            if uv is not None:
                for l, (a, b) in zip(f.loops, ids):
                    l[uvl].uv = uv[a][b]
            if expect is not None and (i, k) == tuple(expect[:2]):
                probe = f
    for cap, ri, ci in ((cap0, 0, -1), (cap1, nr - 1, nr - 1)):
        if cap is None:
            continue
        cv = bm.verts.new(cap)
        mi = mb.mi(ml[0] if ri == 0 else ml[-1])
        for k in range(span):
            if skip is not None and skip(ci, k):
                continue
            k2 = (k + 1) % m
            vs = [V[ri][k], V[ri][k2], cv] if ri == 0 else [V[ri][k2], V[ri][k], cv]
            f = bm.faces.new(vs)
            f.material_index = mi
            if uv is not None:
                for l in f.loops:
                    l[uvl].uv = uv[ri][k]
    loose = [v for v in bm.verts if not v.link_faces]
    if loose:
        bmesh.ops.delete(bm, geom=loose, context='VERTS')
    bm.normal_update()
    if expect is not None:
        assert probe is not None, "orientation probe face was skipped"
        if probe.normal.dot(Vector(expect[2])) < 0:
            bmesh.ops.reverse_faces(bm, faces=bm.faces[:])
            bm.normal_update()
    mb.merge(bm, None)


def loft(mb, secs, mat, xf=None):
    """Beam along +X through chamfered-rectangle sections [(x, half width, z0, z1, chamfer)]: 8 quads per bay,
    3 quads per end (the docking arm's loft without pockets)."""
    bm = bmesh.new()
    rings = []
    for (x, w, z0, z1, c) in secs:
        pts = [(-w + c, z0), (w - c, z0), (w, z0 + c), (w, z1 - c), (w - c, z1), (-w + c, z1), (-w, z1 - c), (-w, z0 + c)]
        rings.append([bm.verts.new((x, y, z)) for y, z in pts])
    for i in range(len(rings) - 1):
        for k in range(8):
            k2 = (k + 1) % 8
            bm.faces.new((rings[i][k], rings[i + 1][k], rings[i + 1][k2], rings[i][k2]))
    for r in (rings[0], rings[-1]):
        bm.faces.new((r[7], r[0], r[1], r[2])); bm.faces.new((r[7], r[2], r[3], r[6])); bm.faces.new((r[6], r[3], r[4], r[5]))
    bm.normal_update()
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    mb.merge(bm, mat, xf)


def foil(mb, sections, mat, xf=None):
    """A thin lofted plate (fin, blade, vane) through sections [(centre Vector, chord dir Vector, normal dir Vector,
    chord, thickness)]: an 8-point lens section, quads along, 3 quads on each end."""
    bm = bmesh.new()
    rings = []
    prof = [(-0.50, 0.0), (-0.40, 0.55), (-0.10, 1.0), (0.25, 0.80), (0.50, 0.0), (0.25, -0.80), (-0.10, -1.0), (-0.40, -0.55)]
    for (c, cd, nd, ch, th) in sections:
        rings.append([bm.verts.new(c + cd * (u * ch) + nd * (v * th * 0.5)) for u, v in prof])
    for i in range(len(rings) - 1):
        for k in range(8):
            k2 = (k + 1) % 8
            bm.faces.new((rings[i][k], rings[i + 1][k], rings[i + 1][k2], rings[i][k2]))
    for r in (rings[0], rings[-1]):
        bm.faces.new((r[0], r[1], r[2], r[7])); bm.faces.new((r[7], r[2], r[3], r[6])); bm.faces.new((r[6], r[3], r[4], r[5]))
    bm.normal_update()
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    mb.merge(bm, mat, xf)


class Surface:
    """BVH of a closed mesh with smooth (interpolated) normals: ray hits and nearest points."""
    def __init__(self, me):
        self.co = [v.co.copy() for v in me.vertices]
        self.vn = [v.normal.copy() for v in me.vertices]
        self.polys = [tuple(p.vertices) for p in me.polygons]
        self.bvh = BVHTree.FromPolygons(self.co, self.polys)

    def _n(self, p, fi):
        ids = self.polys[fi]
        w = poly_3d_calc([self.co[i] for i in ids], p)
        n = Vector()
        for wi, i in zip(w, ids):
            n += self.vn[i] * wi
        return n.normalized()

    def ray(self, o, d):
        hit, nrm, fi, dist = self.bvh.ray_cast(Vector(o), Vector(d).normalized())
        if hit is None:
            return None, None
        return hit, self._n(hit, fi)

    def near(self, p):
        hit, nrm, fi, dist = self.bvh.find_nearest(Vector(p))
        return hit, self._n(hit, fi)


def prune_inside(mb, inside_depth, keep_depth=0.015):
    """drop faces of an attached part whose every vertex lies deeper than keep_depth inside the hull"""
    bm = mb.bm
    dead = [f for f in bm.faces if all(inside_depth(v.co) > keep_depth for v in f.verts)]
    bmesh.ops.delete(bm, geom=dead, context='FACES')
    loose = [v for v in bm.verts if not v.link_faces]
    if loose:
        bmesh.ops.delete(bm, geom=loose, context='VERTS')
    return len(dead)


def resample_profile(prof, mats, max_step):
    """split every profile segment longer than max_step into equal pieces (materials follow their segment)"""
    ml = list(mats) if isinstance(mats, (list, tuple)) else [mats] * (len(prof) - 1)
    out, om = [prof[0]], []
    for i in range(len(prof) - 1):
        (r0, a0), (r1, a1) = prof[i], prof[i + 1]
        n = max(1, int(math.ceil(math.hypot(r1 - r0, a1 - a0) / max_step)))
        for j in range(1, n + 1):
            out.append((r0 + (r1 - r0) * j / n, a0 + (a1 - a0) * j / n))
            om.append(ml[i])
    return out, om


def bm_tree(bm):
    bm.verts.ensure_lookup_table()
    return BVHTree.FromBMesh(bm)


def overlaps(bm_a, bm_b):
    """intersecting face pairs between two bmeshes (indices into their face sequences)"""
    bm_a.faces.ensure_lookup_table(); bm_b.faces.ensure_lookup_table()
    return BVHTree.FromBMesh(bm_a).overlap(BVHTree.FromBMesh(bm_b))

# ====================================================================== from mk_topo.py
# Mako (copied from the Spearfish): clean-topology helpers (exec'd after lib2.py and tender_lib.py). Everything here makes quads only:
# caps are Coons-patch grids instead of fans or n-gons, spheres are cube spheres, patches are ladders.
import bmesh, math
from mathutils import Vector, Matrix, Euler


def coons_cap(bm, ring, mat_index=0):
    """fill a closed ring of 4n BMVerts with an n x n quad grid (a Coons patch over its four sides); returns the faces.
    The ring's order sets the faces' winding (same direction as the ring)."""
    m = len(ring)
    assert m % 4 == 0 and m >= 4, m
    n = m // 4
    P = [v.co.copy() for v in ring]
    B = [P[i] for i in range(0, n + 1)]                      # u 0..1 along side 0
    R = [P[n + i] for i in range(0, n + 1)]                  # v 0..1 along side 1
    T = [P[(3 * n - i) % m] for i in range(0, n + 1)]        # u 0..1, side 2 reversed
    L = [P[(4 * n - i) % m] for i in range(0, n + 1)]        # v 0..1, side 3 reversed
    G = [[None] * (n + 1) for _ in range(n + 1)]             # G[i][j], i along u, j along v
    for i in range(n + 1):
        G[i][0] = ring[i]
        G[i][n] = ring[(3 * n - i) % m]
    for j in range(n + 1):
        G[0][j] = ring[(4 * n - j) % m]
        G[n][j] = ring[n + j]
    for i in range(1, n):
        for j in range(1, n):
            u, v = i / n, j / n
            p = ((1 - v) * B[i] + v * T[i] + (1 - u) * L[j] + u * R[j]
                 - ((1 - u) * (1 - v) * B[0] + u * (1 - v) * B[n] + (1 - u) * v * T[0] + u * v * T[n]))
            G[i][j] = bm.verts.new(p)
    faces = []
    for i in range(n):
        for j in range(n):
            f = bm.faces.new((G[i][j], G[i + 1][j], G[i + 1][j + 1], G[i][j + 1]))
            f.material_index = mat_index
            faces.append(f)
    return faces


def ring_rows(bm, rows, closed=True, mat_index=0):
    """quads between consecutive rows of BMVerts (equal length)"""
    out = []
    for a, b in zip(rows, rows[1:]):
        n = len(a)
        for k in range(n if closed else n - 1):
            k2 = (k + 1) % n
            f = bm.faces.new((a[k], b[k], b[k2], a[k2]))
            f.material_index = mat_index
            out.append(f)
    return out


def orient_out(bm, centre):
    """turn every face so it points away from centre (closed, roughly convex parts)"""
    bm.normal_update()
    flip = [f for f in bm.faces if f.normal.dot(f.calc_center_median() - Vector(centre)) < 0]
    if flip:
        bmesh.ops.reverse_faces(bm, faces=flip)
    bm.normal_update()


def qsphere(mb, ctr, radii, mat, cuts=3, rot=(0, 0, 0)):
    """a cube sphere (quads only), scaled to radii"""
    bm = bmesh.new()
    bmesh.ops.create_cube(bm, size=2.0)
    bmesh.ops.subdivide_edges(bm, edges=bm.edges[:], cuts=cuts, use_grid_fill=True)
    for v in bm.verts:
        d = v.co.normalized()
        v.co = Vector((d.x * radii[0], d.y * radii[1], d.z * radii[2]))
    bm.normal_update()
    mb.merge(bm, mat, Matrix.Translation(ctr) @ Matrix.LocRotScale(None, Euler(rot), None))


def revolve(mb, prof, mats, segs=16, xf=None, cap0=True, cap1=True, cap0_mat=None, cap1_mat=None):
    """revolve a profile [(r, a)] (r > 0 everywhere) about local +Z; open ends closed by Coons caps (segs % 4 == 0)"""
    bm = bmesh.new()
    rows = [[bm.verts.new((r * math.cos(2 * math.pi * k / segs), r * math.sin(2 * math.pi * k / segs), a)) for k in range(segs)]
            for r, a in prof]
    ml = list(mats) if isinstance(mats, (list, tuple)) else [mats] * (len(prof) - 1)
    for i in range(len(rows) - 1):
        mi = mb.mi(ml[i])
        ring_rows(bm, [rows[i], rows[i + 1]], mat_index=mi)
    if cap0:
        coons_cap(bm, rows[0][::-1], mb.mi(cap0_mat or ml[0]))
    if cap1:
        coons_cap(bm, rows[-1], mb.mi(cap1_mat or ml[-1]))
    c = Vector((0, 0, 0.5 * (prof[0][1] + prof[-1][1])))
    bm.normal_update()
    # orientation: the side faces point away from the axis
    flip = []
    for f in bm.faces:
        p = f.calc_center_median()
        radial = Vector((p.x, p.y, 0))
        ref = radial if radial.length > 1e-6 and abs(f.normal.z) < 0.9 else Vector((0, 0, p.z - c.z))
        if f.normal.dot(ref) < 0:
            flip.append(f)
    if flip:
        bmesh.ops.reverse_faces(bm, faces=flip)
    bm.normal_update()
    mb.merge(bm, None, xf)


def qcyl(mb, r, h, ctr, mat, rot=(0, 0, 0), segs=12, r2=None):
    """a cylinder along local Z with Coons-cap ends (quads only)"""
    r2 = r if r2 is None else r2
    xf = Matrix.Translation(ctr) @ Matrix.LocRotScale(None, Euler(rot), None)
    revolve(mb, [(r, -h / 2), (r2, h / 2)], mat, segs=segs, xf=xf)


def qbox(mb, size, ctr, mat, rot=(0, 0, 0), bevel=0.0):
    """a box; bevelled with two segments, which keeps the corners quads"""
    mb.box(size, ctr, mat, rot=rot, bevel=bevel, seg=2 if bevel > 0 else 1)


def ladder_patch(mb, lo_hi, t0, t1, to3d, mat, nt=8, nacross=3, lift=0.002, bulge=0.0, normal=None):
    """a thin plate over a region given as t in [t0, t1] -> (s_lo, s_hi) across: a grid of nt x nacross quads.
    to3d(t, s) -> point on the hull; normal(p) -> outward normal. Ends at t0/t1 keep their width (no fan)."""
    bm = bmesh.new()
    rows = []
    for i in range(nt + 1):
        t = t0 + (t1 - t0) * i / nt
        lo, hi = lo_hi(t)
        row = []
        for j in range(nacross + 1):
            s = lo + (hi - lo) * j / nacross
            p = to3d(t, s)
            n = normal(p)
            w = 1.0 - (2.0 * j / nacross - 1.0) ** 2
            e = 1.0 - (2.0 * i / nt - 1.0) ** 2
            row.append(bm.verts.new(p + n * (lift + bulge * max(0.0, w) * max(0.0, e) ** 0.5)))
        rows.append(row)
    ring_rows(bm, rows, closed=False)
    bm.normal_update()
    mid = rows[nt // 2][nacross // 2].co
    nm = normal(mid)
    if sum((f.normal.dot(nm) for f in bm.faces)) < 0:
        bmesh.ops.reverse_faces(bm, faces=bm.faces[:])
    mb.merge(bm, mat)
    return rows


def stadium(t0, t1, half):
    """lo_hi for a stadium (capsule) from t0 to t1 of half width `half`: round ends"""
    r = half
    def f(t):
        a, b = min(t0, t1), max(t0, t1)
        if t < a + r:
            d = (a + r - t) / r
        elif t > b - r:
            d = (t - (b - r)) / r
        else:
            d = 0.0
        w = r * math.sqrt(max(0.0, 1.0 - min(1.0, d) ** 2))
        return -max(w, 0.18 * r), max(w, 0.18 * r)
    return f