#! /usr/bin/python3
# Last edited on 2024-07-27 00:30:06 by stolfi

from math import sin, cos, tan, atan2, asin, hypot, sqrt, pi, inf, floor, sqrt
import rn, rmxn
import sys
import re
import random

def build():

  Nrack = 20
  Npini = 20
  wd = 20  # Width of rack and pinion.
  th = 10  # Pitch (teeth spacing along datum line)
  
  rack = build_rack(Nrack,th,wd,bw)
  pinion = build_pinion(Npini,th,wd,bw)
  
  V, E, F = mesh_from_vertices(vtot_rack, vtot_pini);
  M = { 
    'V': rack['V'] + pinion['V'],
    'E': rack['E'] + pinion['E'],
    'F': rack['F'] + pinion['F'],
    'T': rack['T'] + pinion['T'],
  }
  return M
  
def build_rack(nt,nf,th):
  # Returns a Python dict {M} with the mesh of the rack.
  #
  # The rack will have {nt} teeth with pitch {th} along the datum line.
  # The long direction is {Z}, and the teeth are on the {+X} face. It is
  # basically a prism obtained by taking profile polygon on the {Y=0} plane,
  # extruding it in {Y} by {±wd/2}, then shaving off a 45 degree bevel
  # with depth {bw} along the border of each base. 
  # 
  # Each tooth will have the generic tooth profile as defined by
  # {build_tooth_profile(nf,0,nf,th)}.
  #
  P = rack_profile(nt,nf,th)
  M = extrude_profile(P,wd,bw)
  M['T'] = rack_triang_edges(M['V'])
  return M

def rack_profile(nt,nf,th):
  # Returns a list {vp} with the vertices of the rack's profile polygon. 
  #
  # The rack will have {nt} teeth with pitch {th} along the datum line.
  # The profile polygon lies on the plane {Y=0}. The long direction is
  # {Z}, and the teeth are pointing in the {+X} side.
  # 
  # Each tooth will have the generic tooth profile as defined by
  # {build_tooth_profile(nf,0,nf,th)}.
  #
  # Each vertex is a 4-tuple where the first three elements are the 
  # coordinates {x,y,z} (with {y=0}) and the fourth element is a string label.
  # The label is "B.{KKKK}" for the straight part of the
  # profile and "T.{IIII}.{u}.{JJJJ}" for the vertices of the teeth.
  #
  # The index {k} for the 'B' vertices (the corners of the flat face)
  # can be {0..5}. The vertices with {k=0} and {k=5} coincide with the
  # bottom vertices of the first and last teeth; see below.
  #
  # The index {i} for the 'T' vertices is the tooth index, in {0..nt-1},
  # from low {Z} to high {Z}. The key {u} is 'b0', 't0', 't1', 'b1'
  # where the first letter is 'b' for the bottom of the tooth gap, 't'
  # for the tooth tip, and the second letter is '0' for the low-{Z} half
  # of the tooth, '1' for the high-{Z} half. The index {j} ranges in
  # {0..nf}, in increasing {Z} for {u} = 'b0' or 't1' and decreasing {Z}
  # for {u} = 't0' or 'b1'.
  #
  # Some of the vertices coincide and only one of each coincident set
  # should be used in the mesh. Specifically, for all {s} and all {i} in
  # {0..nt-1}:
  #
  #   {vtot['B'][s][0]} = {vtot['T'][s][0]['b0']}.
  #
  #   {vtot['B'][s][5]} = {vtot['T'][s][nt-1]['b1']}.
  #
  #   {vtot['T'][s][i]['t0'][0]} = {vtot['T'][s][i]['t1'][0]}
  #   
  #   {vtot['T'][s][i]['b1'][0]} = {vtot['T'][s][i+1]['b0'][0]}
  # like {s}
  # is 'pf', 'pe', 'me', or 'mf', where the first letter specifies the
  # sign of the {Y} coordinate (positive 'p', negative 'm'), and the
  # second is 'f' for the flat face at {Y = ±wd/2} and 'e' for the bevel
  # edge at {Y = ±(wd/2 - bw)}.

def extrude_and_bevel(vref,wd,bw):
  # Returns a Python dict {vtot} with all the vertices of the rack. 
  #
  # The rack will have {nt} teeth with pitch {th} along the datum line.
  # The long direction is {Z}, and the teeth are on the {+X} face. It is
  # basically a prism obtained by taking a /profile/ polygon on the
  # {Y=0} plane, extruding it in {Y} by {±wd/2}, then shaving off a 45
  # degree bevel with depth {bw} along the border of each base. Each
  # tooth will have {4*nf+2} rectangular faces paralle to the {Y} axis,
  # not counting the bevel faces.
  #
  # The vertices are {vtot[s][k]} where {k} varies from {0..Nu-1} where
  # {Nu} is the number of vertices in the profile polygon. The key {s}
  # is 'pf', 'pe', 'me', or 'mf', where the first letter specifies the
  # sign of the {Y} coordinate (positive 'p', negative 'm'), and the
  # second is 'f' for the flat face at {Y = ±wd/2} and 'e' for the bevel
  # edge at {Y = ±(wd/2 - bw)}.
  # 
  # For both sets, the key {s} is 'pf', 'pe', 'me', or 'mf', where the
  # first letter specifies the sign of the {Y} coordinate (positive 'p',
  # negative 'm'), and the second is 'f' for the flat face at {Y =
  # ±wd/2} and 'e' for the bevel edge at {Y = ±(wd/2 - bw)}.

  vtot = { 'mf': [], 'me': [], 'pe': [], 'pf': [] }
  
  # Basic profile polygon, before extrusion and beveling:
  vref = rack_vertices_ref_dict(nt)
  
  Nu = len(vref)
  
  for k in range(Nu):
    k_prev = (k - 1) % Nu
    k_next = (k + 1) % Nu
    x_prev, z_prev = vref[k_prev][0], vref[k_prev][2]
    x_this, z_this = vref[k][0], vref[k][2]
    x_next, z_next = vref[k_next][0], vref[k_next][2]
    for dy in 0, 1:
      sy = "mp"[dy]
      yf = 0.5*wd*(2*dy - 1)
      ye = (0.5*wd - bw)*(2*dy - 1)
      vtot[sy+'e'].append ( x_this,ye,z_this, 'R.'+sy+'e'.vref[k][3] )
      
      # Compute the bevel vertex on the face:
      x_beve,z_beve = compute_bevel_vertex(x_prev,z_prev, x_this,z_this, x_next,z_next, bw)
      vtot['B'][sy+'f'][k] = ( x_beve,yf,z_beve, 'R.'+sy+'f'.vref[k][3] )
      
        s = "mp"[sy] + "fe"[sb]
        for sz in +1, 0, -1:
          key_src = f"{'pop'[sx+1]}{'pop'[sy+1]}{'pop'[sz+1]}"
          v_src = vref.get(key_src, None)
          if v_src != None and len(v_src) != 0:
            Nv_src = len(v_src)
            key_dst = f"{'mop'[sx+1]}{'mop'[sy+1]}{'mop'[sz+1]}"
            v_dst = [ None ] * Nv_src
            for iv in range(Nv_src):
              Nv += 1
              vi_src = v_src[iv]
              vi_dst = ( vi_src[0]*sx, vi_src[1]*sy, vi_src[2]*sz, Nv, ) 
              v_dst[iv] = vi_dst
            vtot[key_dst] = v_dst
  return vtot, Nv

def ref_vertices(S):
  # Computes the vertices of one eight of the solid. Parameters:
  
  # {S}        Side of main cube.
  
  # The value of {S} is assumed approximate. The actual dimensions
  # will be rounded to apropriate even multiples of {eps}.

  # Returns a dict with keys 'vopp', 'vpop', 'vppo', 'vppp'.
  # The values are arrays of vertices.
  # Array 'vppp' is the vertices in the positive octant.
  # Array 'vopp' has the vertices on positive quadrant of {YZ} plane.  
  # Array 'vpop' has the vertices on positive quadrant of {XZ} plane.  
  # Array 'vppo' has the vertices on positive quadrant of {XY} plane.  
  #
  # Each element is a triplet of coordinates. 
  # Does not round any coordinate.

  sys.stderr.write("ref_vertices:")
  sys.stderr.write(" S = %.*f \n" % (prec, S))
  
  # Approximate key dimensions:
  HS = S/2        # Half-side of main cube.
  Ds = HS/4       # Depth of main slot in main cube.
  Es = HS/4       # Distance in {X} from slot to edge of main cube.
  Lh = HS/2       # Length in {X} of handle.
  Wh = HS/3       # Half-width in {Y} of handle.
  Th = HS/4       # Half-height in {Z} of handle in 
  Eh = Wh/2       # Thickness of handle all around. 
  
  Xk = [ None ] * 24              
  Xk[ 0] = 1.0000 * HS      # {X} of side of main cube.
  Xk[ 1] = 0.2500 * HS      # Low {X} of slot and holes on main cube.
  Xk[ 2] = 0.7500 * HS      # High {X} of slot and holes on main cube.
  Xk[ 3] = 1.6000 * HS      # High {X} of side handle.
  Xk[ 4] = 1.2000 * HS      # Low {X} of hole in handle.
  Xk[ 5] = 1.4000 * HS      # High {X} of hole in handle. 
  Xk[ 6] = 0.4000 * HS      # Low {X} of bar in hole.
  Xk[ 7] = 0.6000 * HS      # High {X} of bar in hole.
  Xk[ 8] = 0.6250 * HS      # Some {X} of pockmarks.
  Xk[ 9] = 0.6500 * HS      # Some {X} of pockmarks.
  Xk[10] = 0.6750 * HS      # Some {X} of pockmarks.
  Xk[11] = 0.7000 * HS      # Some {X} of pockmarks.
  Xk[12] = 0.7250 * HS      # Some {X} of pockmarks.
  Xk[13] = 0.7750 * HS      # Some {X} of pockmarks.
  Xk[14] = 0.8000 * HS      # Some {X} of pockmarks.
  Xk[15] = 0.8250 * HS      # Some {X} of pockmarks.
  Xk[16] = 0.8500 * HS      # Some {X} of pockmarks.
  Xk[17] = 0.8750 * HS      # Some {X} of pockmarks.
  Xk[18] = 0.9000 * HS      # Some {X} of pockmarks.
  Xk[19] = 0.9250 * HS      # Some {X} of pockmarks.
  Xk[20] = 0.6625 * HS      # Some {X} of pockmarks.
  Xk[21] = 0.7125 * HS      # Some {X} of pockmarks.
  Xk[22] = 0.8375 * HS      # Some {X} of pockmarks.
  Xk[23] = 0.8875 * HS      # Some {X} of pockmarks.
 
  Yk = [ None ] * 4
  Yk[ 0] = 1.0000 * HS      # {Y} of side  of main cube.
  Yk[ 1] = 0.3000 * HS      # High {Y} of handle.
  Yk[ 2] = 0.1000 * HS      # High {Y} of hole in handle.
  Yk[ 3] = 0.9500 * HS      # {Y} of bottom of pockmarks.
  
  Zk = [ None ] * 22
  Zk[ 0] = 1.0000 * HS      # High {Z} of main cube.
  Zk[ 1] = 0.8000 * HS      # {Z} of bottom of slot on main cube.
  Zk[ 2] = 0.5500 * HS      # {Z} of top of hole on main cube.
  Zk[ 3] = 0.1500 * HS      # {Z} of top of handle.
  Zk[ 4] = 0.1500 * HS      # {Z} of bottom of bar in hole.
  Zk[ 5] = 0.4000 * HS      # {Z} of top of bar in hole.
  Zk[ 6] = 0.4250 * HS      # Some {Z} of pockmarks.
  Zk[ 7] = 0.4500 * HS      # Some {Z} of pockmarks.
  Zk[ 8] = 0.4750 * HS      # Some {Z} of pockmarks.
  Zk[ 9] = 0.5000 * HS      # Some {Z} of pockmarks.
  Zk[10] = 0.5250 * HS      # Some {Z} of pockmarks.
  Zk[11] = 0.5750 * HS      # Some {Z} of pockmarks.
  Zk[12] = 0.6000 * HS      # Some {Z} of pockmarks.
  Zk[13] = 0.6250 * HS      # Some {Z} of pockmarks.
  Zk[14] = 0.6500 * HS      # Some {Z} of pockmarks.
  Zk[15] = 0.6750 * HS      # Some {Z} of pockmarks.
  Zk[16] = 0.7000 * HS      # Some {Z} of pockmarks.
  Zk[17] = 0.7250 * HS      # Some {Z} of pockmarks.
  Zk[18] = 0.4625 * HS      # Some {Z} of pockmarks.
  Zk[19] = 0.5125 * HS      # Some {Z} of pockmarks.
  Zk[20] = 0.6375 * HS      # Some {Z} of pockmarks.
  Zk[21] = 0.6875 * HS      # Some {Z} of pockmarks.
  
  vopp = [ None ] *  2
  vpop = [ None ] *  0
  vppo = [ None ] *  0
  vppp = [ None ] * (15 + (28 if pock else 0))
  
  vopp[ 0] = (0,     Yk[ 0], Zk[ 0])
  vopp[ 1] = (0,     Yk[ 0], Zk[ 2])

  vppp[ 0] = (Xk[ 1], Yk[ 0], Zk[ 0])
  vppp[ 1] = (Xk[ 1], Yk[ 0], Zk[ 1])
  vppp[ 2] = (Xk[ 2], Yk[ 0], Zk[ 1])
  vppp[ 3] = (Xk[ 2], Yk[ 0], Zk[ 0])
  vppp[ 4] = (Xk[ 0], Yk[ 0], Zk[ 0])
  vppp[ 5] = (Xk[ 1], Yk[ 0], Zk[ 2])
  vppp[ 6] = (Xk[ 2], Yk[ 0], Zk[ 2])
  vppp[ 7] = (Xk[ 0], Yk[ 1], Zk[ 3])
  vppp[ 8] = (Xk[ 3], Yk[ 1], Zk[ 3])
  vppp[ 9] = (Xk[ 4], Yk[ 2], Zk[ 3])
  vppp[10] = (Xk[ 5], Yk[ 2], Zk[ 3])

  #  Bars in holes
  vppp[11] = (Xk[ 6], Yk[ 0], Zk[ 5])
  vppp[12] = (Xk[ 7], Yk[ 0], Zk[ 5])
  vppp[13] = (Xk[ 6], Yk[ 0], Zk[ 4])
  vppp[14] = (Xk[ 7], Yk[ 0], Zk[ 4])
  
  if pock:
    # Pockmark top corners:
    vppp[15] = (Xk[ 2], Yk[ 0], Zk[12])
    vppp[16] = (Xk[11], Yk[ 0], Zk[12])
    vppp[17] = (Xk[11], Yk[ 0], Zk[14])
    vppp[18] = (Xk[10], Yk[ 0], Zk[12])
    vppp[19] = (Xk[ 8], Yk[ 0], Zk[13])
    vppp[20] = (Xk[ 9], Yk[ 0], Zk[15])
    vppp[21] = (Xk[10], Yk[ 0], Zk[16])
    vppp[22] = (Xk[12], Yk[ 0], Zk[17])
    vppp[23] = (Xk[ 2], Yk[ 0], Zk[15])
    vppp[24] = (Xk[ 2], Yk[ 0], Zk[14])
    vppp[25] = (Xk[14], Yk[ 0], Zk[12])
    vppp[26] = (Xk[14], Yk[ 0], Zk[ 2])
    vppp[27] = (Xk[16], Yk[ 0], Zk[ 2])
    vppp[28] = (Xk[16], Yk[ 0], Zk[ 9])
    vppp[29] = (Xk[17], Yk[ 0], Zk[ 2])
    vppp[30] = (Xk[19], Yk[ 0], Zk[10])
    vppp[31] = (Xk[18], Yk[ 0], Zk[ 8])
    vppp[32] = (Xk[17], Yk[ 0], Zk[ 7])
    vppp[33] = (Xk[15], Yk[ 0], Zk[ 6])
    vppp[34] = (Xk[14], Yk[ 0], Zk[ 8])
    vppp[35] = (Xk[14], Yk[ 0], Zk[ 9])

    # Pockmark bottoms:
    vppp[36] = (Xk[15], Yk[ 3], Zk[10])
    vppp[37] = (Xk[12], Yk[ 3], Zk[13])
    vppp[38] = (Xk[22], Yk[ 3], Zk[18])
    vppp[39] = (Xk[23], Yk[ 3], Zk[19])
    vppp[40] = (Xk[20], Yk[ 3], Zk[20])
    vppp[41] = (Xk[21], Yk[ 3], Zk[21])
    vppp[42] = (Xk[13], Yk[ 3], Zk[11])

  return { 'opp': vopp, 'pop': vpop, 'ppo': vppo, 'ppp': vppp }
  
def make_faces_OBJ(wro, vtot, pock):
  # Creates the faces of the object to OBJ file {wro}, assuming that
  # the vertices are in the dict {vtot}.
  
  # Also returns tables {Fnrm} of normals to the faces and {Fbar} of their
  # barycenters.
  #
  # Also returns tables {E,G} of edges:
  #
  #   {E} Non-ghost edges of the mesh.
  #   {G} Ghost edges of the mesh.
  #
  # Each element of these arrays is a pair {(iorg,idst)} of OBJ indices
  # of the endpoint vertices. Each array {Fnrm,Fbar,E,G} is indexed
  # {[1..N]}; element 0 is not used ({None}).
  
  # Expected face, and edge counts perpendicular or parallel to each axis, and
  # general oblique:
  
  #        main shape      
  #        -----------------
  #           X     Y     Z
  #        -----------------       
  Nf_exp =   28 +  18 +  26   # Expected face counts (perp).
  Ne_exp =   64 +  60 +  60   # Expected non-ghost edge count (para).
  Ng_exp =    2 +   6 +   4   # Expected ghost edge count (perp).
  
  if pock:
    #           ----------------------
    #              X     Y     Z  oblq
    #           ----------------------
    Nf_exp +=      0 +   0 +   0 + 224
    Ne_exp +=     48 +   0 +  48 + 288
    Ng_exp +=      0 +   0 +   0 +   0
    
  # Variables that are global to the nested functions:

  E = [ None ]
  G = [ None ]
  Fbar = [ None ]
  Fnrm = [ None ]

  debug_face = False
  Fp = []    # Vertices of current face, indexed from 0.
  
  def Sved(p_org, p_dst, ghost):
    # Saves the edge from {p_org} to {p_dst} in {E} or {G} depending on
    # {ghost}, if traversed in increasing vertex index sense. 

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    kv_org = p_org[3]
    kv_dst = p_dst[3]
    # Avoid duplicate edges (even ghost ones):
    if kv_org < kv_dst:
      if ghost:
        G.append((kv_org, kv_dst))
      else:
        E.append((kv_org, kv_dst))
    return None

  def Ptit(tit):
    # Title for a face or set of faces.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    wro.write("\n");
    wro.write("# %s\n" % tit);
    sys.stderr.write("%s \n" % (tit));
    return None
  
  def Bof(debug = False):
    #  Start of a new face.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    debug_face = debug
    assert len(Fp) == 0
    return None
  
  def Prtv(key, iv, ghost):
    # Adds vertex {vtot[key][iv]} to the current face, augmented with 
    # its "{key}:{iv}" label.  
    #
    # If {ghost} is {True}, assumes that the edge from previous
    # vertex is a ghost one.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face
    p = vtot[key][iv]
    assert len(p) == 4
    # Saves vertex in list {Fp} of face vertices, with label:
    lab = f"{key}:{iv:02}"
    Fp.append(p + ( lab, ));
    if len(Fp) >= 2: Sved(Fp[-2], Fp[-1], ghost)
    return None
    
  def Eof(ghost, order):
    # End of face. If {ghost}, assumes that 
    # closing edge is a ghost one.
    # Prints in given or reverse order depending on {order} is {+1} or {-1}.

    nonlocal E, G, Fbar, Fnrm, Fp, debug_face

    if debug_face:
      sys.stderr.write("  writing face %d (%d corners)\n" % (len(Fnrm), len(Fp))) 
    assert len(Fp) >= 3
    Sved(Fp[-1], Fp[0], ghost)

    # Write the face:
    wro.write("f")
    deg = len(Fp)
    for rv in range(deg):
      jv = rv if order == +1 else deg-1-rv
      p = Fp[jv]
      assert len(p) == 5
      kv = p[3]
      lab = p[4]
      if debug_face:
        sys.stderr.write("    v%04d = %s = ( %9.*f %9.*f %9.*f )\n" % (kv, lab, prec, p[0], prec, p[1], prec, p[2]));
      wro.write(" %d" % kv);
    wro.write("\n");
    
    # Compute normal, barycenter:
    nrm, bar = face_normal_and_barycenter(Fp)
    Fnrm.append(nrm)
    Fbar.append(bar)
    Fp = []; 
    return None
    
  # Here are the faces:
  
  debug_faces = True

  Ptit(f"faces with normal ±X")
  for sx in +1, -1:
    tx = 'mop'[sx+1]

    Ptit(f"  faces in halfspace sx = {sx:+d}")
    
    # Vert face of cube:
    Bof(debug_faces)
    Prtv(tx+'pp',  4, False)
    
    Prtv(tx+'mp',  4, False)
    Prtv(tx+'mm',  4, False)
    Prtv(tx+'pm',  4, False)
    Prtv(tx+'pp',  4, False)
    
    Prtv(tx+'pp',  7, True)
    Prtv(tx+'pm',  7, False)
    Prtv(tx+'mm',  7, False)
    Prtv(tx+'mp',  7, False)
    Prtv(tx+'pp',  7, False)
    Eof(True, sx)
    
    # Centrad vert wall of hole on cube:
    Bof(debug_faces)
    Prtv(tx+'pp',  5, False)
    Prtv(tx+'mp',  5, False)
    Prtv(tx+'mm',  5, False)
    Prtv(tx+'pm',  5, False)
    Eof(False, sx)
    
    # Distad vert wall of hole on cube:
    Bof(debug_faces)
    Prtv(tx+'pp',  6, False)
    Prtv(tx+'pm',  6, False)
    Prtv(tx+'mm',  6, False)
    Prtv(tx+'mp',  6, False)
    Eof(False, sx)
    
    # Centrad vert wall of hole on handle:
    Bof(debug_faces)
    Prtv(tx+'pp',  9, False)
    Prtv(tx+'mp',  9, False)
    Prtv(tx+'mm',  9, False)
    Prtv(tx+'pm',  9, False)
    Eof(False, sx)
    
    # Distad vert wall of hole on handle:
    Bof(debug_faces)
    Prtv(tx+'pp', 10, False)
    Prtv(tx+'pm', 10, False)
    Prtv(tx+'mm', 10, False)
    Prtv(tx+'mp', 10, False)
    Eof(False, sx)
    
    # Extreme vert wall of handle:
    Bof(debug_faces)
    Prtv(tx+'pp',  8, False)
    Prtv(tx+'mp',  8, False) 
    Prtv(tx+'mm',  8, False)
    Prtv(tx+'pm',  8, False)
    Eof(False, sx)

    for sz in +1, -1:
      tz = 'mop'[sz+1]

      Ptit(f"    faces in quadrant sx = {sx:+d}, sz = {sz:+d}")
      
      # Centrad vert wall of slot on face of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  0, False)
      Prtv(tx+'m'+tz,  0, False)
      Prtv(tx+'m'+tz,  1, False)
      Prtv(tx+'p'+tz,  1, False)
      Eof(False, sx*sz)
      
      # Distad vert wall of slot on {±Z} face of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  2, False)
      Prtv(tx+'m'+tz,  2, False)
      Prtv(tx+'m'+tz,  3, False)
      Prtv(tx+'p'+tz,  3, False)
      Eof(False, sx*sz)
       
      # Distad vert wall of bar in hole of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz, 12, False)
      Prtv(tx+'m'+tz, 12, False)
      Prtv(tx+'m'+tz, 14, False)
      Prtv(tx+'p'+tz, 14, False)
      Eof(False, sx*sz)
      
      # Centrad vert wall of bar in hole of cube:
      Bof(debug_faces)
      Prtv(tx+'m'+tz, 11, False)
      Prtv(tx+'p'+tz, 11, False)
      Prtv(tx+'p'+tz, 13, False)
      Prtv(tx+'m'+tz, 13, False)
      Eof(False, sx*sz)

  Ptit(f"faces with normal ±Y")

  for sy in +1, -1:
    ty = 'mop'[sy+1]
    
    Ptit(f"  faces in halfpspace sy = {sy:+d}")
    
    # Vert face of main cube:
    Ptit(f"    big complicated face with holes and pinch corners:")
    
    Prtv('o'+ty+'m',  0, False)
    Prtv('o'+ty+'m',  1, True)
    
    Prtv('p'+ty+'m',  5, True)
    Prtv('p'+ty+'m',  6, False)

    if pock:
      Prtv('p'+ty+'m', 15, False)
      Prtv('p'+ty+'m', 16, False)
      Prtv('p'+ty+'m', 17, False)
      Prtv('p'+ty+'m', 18, False)
      Prtv('p'+ty+'m', 19, False)
      Prtv('p'+ty+'m', 20, False)
      Prtv('p'+ty+'m', 17, False)
      Prtv('p'+ty+'m', 21, False)
      Prtv('p'+ty+'m', 22, False)
      Prtv('p'+ty+'m', 23, False)
      Prtv('p'+ty+'m', 17, False)
      Prtv('p'+ty+'m', 24, False)
      Prtv('p'+ty+'m', 15, False)
      Prtv('p'+ty+'m', 25, False)
      Prtv('p'+ty+'m', 26, False)
      Prtv('p'+ty+'m', 27, False)
      Prtv('p'+ty+'m', 28, False)
      Prtv('p'+ty+'m', 29, False)
      Prtv('p'+ty+'m', 30, False)
      Prtv('p'+ty+'m', 31, False)
      Prtv('p'+ty+'m', 28, False)
      Prtv('p'+ty+'m', 32, False)
      Prtv('p'+ty+'m', 33, False)
      Prtv('p'+ty+'m', 34, False)
      Prtv('p'+ty+'m', 28, False)
      Prtv('p'+ty+'m', 35, False)
      Prtv('p'+ty+'m', 26, False)
      Prtv('p'+ty+'m',  6, False)

    Prtv('p'+ty+'p',  6, False)

    if pock:
      Prtv('p'+ty+'p', 26, False)
      Prtv('p'+ty+'p', 35, False)
      Prtv('p'+ty+'p', 28, False)
      Prtv('p'+ty+'p', 34, False)
      Prtv('p'+ty+'p', 33, False)
      Prtv('p'+ty+'p', 32, False)
      Prtv('p'+ty+'p', 28, False)
      Prtv('p'+ty+'p', 31, False)
      Prtv('p'+ty+'p', 30, False)
      Prtv('p'+ty+'p', 29, False)
      Prtv('p'+ty+'p', 28, False)
      Prtv('p'+ty+'p', 27, False)
      Prtv('p'+ty+'p', 26, False)
      Prtv('p'+ty+'p', 25, False)
      Prtv('p'+ty+'p', 15, False)
      Prtv('p'+ty+'p', 24, False)
      Prtv('p'+ty+'p', 17, False)
      Prtv('p'+ty+'p', 23, False)
      Prtv('p'+ty+'p', 22, False)
      Prtv('p'+ty+'p', 21, False)
      Prtv('p'+ty+'p', 17, False)
      Prtv('p'+ty+'p', 20, False)
      Prtv('p'+ty+'p', 19, False)
      Prtv('p'+ty+'p', 18, False)
      Prtv('p'+ty+'p', 17, False)
      Prtv('p'+ty+'p', 16, False)
      Prtv('p'+ty+'p', 15, False)
      Prtv('p'+ty+'p',  6, False)

    Prtv('p'+ty+'p',  5, False)
    Prtv('p'+ty+'m',  5, False)
    
    Prtv('o'+ty+'m',  1, True)
    
    Prtv('m'+ty+'m',  5, True)
    Prtv('m'+ty+'p',  5, False)
    
    Prtv('m'+ty+'p',  6, False)

    if pock:
      Prtv('m'+ty+'p', 15, False)
      Prtv('m'+ty+'p', 16, False)
      Prtv('m'+ty+'p', 17, False)
      Prtv('m'+ty+'p', 18, False)
      Prtv('m'+ty+'p', 19, False)
      Prtv('m'+ty+'p', 20, False)
      Prtv('m'+ty+'p', 17, False)
      Prtv('m'+ty+'p', 21, False)
      Prtv('m'+ty+'p', 22, False)
      Prtv('m'+ty+'p', 23, False)
      Prtv('m'+ty+'p', 17, False)
      Prtv('m'+ty+'p', 24, False)
      Prtv('m'+ty+'p', 15, False)
      Prtv('m'+ty+'p', 25, False)
      Prtv('m'+ty+'p', 26, False)
      Prtv('m'+ty+'p', 27, False)
      Prtv('m'+ty+'p', 28, False)
      Prtv('m'+ty+'p', 29, False)
      Prtv('m'+ty+'p', 30, False)
      Prtv('m'+ty+'p', 31, False)
      Prtv('m'+ty+'p', 28, False)
      Prtv('m'+ty+'p', 32, False)
      Prtv('m'+ty+'p', 33, False)
      Prtv('m'+ty+'p', 34, False)
      Prtv('m'+ty+'p', 28, False)
      Prtv('m'+ty+'p', 35, False)
      Prtv('m'+ty+'p', 26, False)
      Prtv('m'+ty+'p',  6, False)

    Prtv('m'+ty+'m',  6, False)

    if pock:
      Prtv('m'+ty+'m', 26, False)
      Prtv('m'+ty+'m', 35, False)
      Prtv('m'+ty+'m', 28, False)
      Prtv('m'+ty+'m', 34, False)
      Prtv('m'+ty+'m', 33, False)
      Prtv('m'+ty+'m', 32, False)
      Prtv('m'+ty+'m', 28, False)
      Prtv('m'+ty+'m', 31, False)
      Prtv('m'+ty+'m', 30, False)
      Prtv('m'+ty+'m', 29, False)
      Prtv('m'+ty+'m', 28, False)
      Prtv('m'+ty+'m', 27, False)
      Prtv('m'+ty+'m', 26, False)
      Prtv('m'+ty+'m', 25, False)
      Prtv('m'+ty+'m', 15, False)
      Prtv('m'+ty+'m', 24, False)
      Prtv('m'+ty+'m', 17, False)
      Prtv('m'+ty+'m', 23, False)
      Prtv('m'+ty+'m', 22, False)
      Prtv('m'+ty+'m', 21, False)
      Prtv('m'+ty+'m', 17, False)
      Prtv('m'+ty+'m', 20, False)
      Prtv('m'+ty+'m', 19, False)
      Prtv('m'+ty+'m', 18, False)
      Prtv('m'+ty+'m', 17, False)
      Prtv('m'+ty+'m', 16, False)
      Prtv('m'+ty+'m', 15, False)
      Prtv('m'+ty+'m',  6, False)

    Prtv('m'+ty+'m',  5, False)
    
    Prtv('o'+ty+'m',  1, True)
    Prtv('o'+ty+'m',  0, True)
    
    Prtv('m'+ty+'m',  0, False)
    Prtv('m'+ty+'m',  1, False)
    Prtv('m'+ty+'m',  2, False)
    Prtv('m'+ty+'m',  3, False)
    Prtv('m'+ty+'m',  4, False)

    Prtv('m'+ty+'p',  4, False)
    Prtv('m'+ty+'p',  3, False)
    Prtv('m'+ty+'p',  2, False)
    Prtv('m'+ty+'p',  1, False)
    Prtv('m'+ty+'p',  0, False)
    
    Prtv('o'+ty+'p',  0, False)
    
    Prtv('p'+ty+'p',  0, False)
    Prtv('p'+ty+'p',  1, False)
    Prtv('p'+ty+'p',  2, False)
    Prtv('p'+ty+'p',  3, False)
    Prtv('p'+ty+'p',  4, False)

    Prtv('p'+ty+'m',  4, False)
    Prtv('p'+ty+'m',  3, False)
    Prtv('p'+ty+'m',  2, False)
    Prtv('p'+ty+'m',  1, False)
    Prtv('p'+ty+'m',  0, False)

    Eof(False, sy)
     
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sx = {sx:+d} sy = {sy:+d}")
    
      # Vert outer face of handle:
      Bof(debug_faces)
      Prtv(tx+ty+'p',  7, False)
      Prtv(tx+ty+'p',  8, False)
      Prtv(tx+ty+'m',  8, False)
      Prtv(tx+ty+'m',  7, False)
      Eof(False, sx*sy)
    
      # Vert inner face of handle:
      Bof(debug_faces)
      Prtv(tx+ty+'p', 10, False)
      Prtv(tx+ty+'p',  9, False)
      Prtv(tx+ty+'m',  9, False)
      Prtv(tx+ty+'m', 10, False)
      Eof(False, sx*sy)
      
      for sz in +1, -1:
        tz = 'mop'[sz+1]

        Ptit(f"      faces in octant sx = {sx:+d} sy = {sy:+d} sz = {sz:+d}")

        # Tip of bar in hole:
        Bof(debug_faces)
        Prtv(tx+ty+tz, 11, False)
        Prtv(tx+ty+tz, 12, False)
        Prtv(tx+ty+tz, 14, False)
        Prtv(tx+ty+tz, 13, False)
        Eof(False, sx*sy*sz)

        if pock:
          Ptit(f"      pockmark faces in octant sx = {sx:+d} sy = {sy:+d} sz = {sz:+d}")
          # Each tuple below is center vertex and perimeter vertices in CCW order.
          Pockmarks = (
            ( 36, 26, 35, 28, 27, ),
            ( 37, 15, 24, 17, 16, ),
            ( 38, 28, 34, 33, 32, ),
            ( 39, 28, 31, 30, 29, ),
            ( 40, 17, 20, 19, 18, ),
            ( 41, 17, 23, 22, 21, ),
            ( 42,  6, 26, 25, 15, ),
          )
          for Pck in Pockmarks:
            ivt = Pck[0]
            for ka in range(4):
              kb = (ka + 1) % 4
              iva = Pck[ka+1]
              ivb = Pck[kb+1]
              Bof(debug_faces)
              Prtv(tx+ty+tz, ivt, False)
              Prtv(tx+ty+tz, iva, False)
              Prtv(tx+ty+tz, ivb, False)
              Eof(False, sx*sy*sz)

  Ptit(f"faces with normal ±Z")

  for sz in +1, -1:
    tz = 'mop'[sz+1]
    
    Ptit(f"  faces in halfspace sz = {sz:+d}")
    
    # Centrad horz face of main cube:
    Bof(debug_faces)
    Prtv('op'+tz,  0, False)
    Prtv('mp'+tz,  0, False)
    Prtv('mm'+tz,  0, False)
    Prtv('om'+tz,  0, False)
    Prtv('pm'+tz,  0, False)
    Prtv('pp'+tz,  0, False)
    Eof(False, sz)
    
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sz = {sz:+d}, sx = {sx:+d}")

      # Distad horz face of main cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  3, False)
      Prtv(tx+'m'+tz,  3, False)
      Prtv(tx+'m'+tz,  4, False)
      Prtv(tx+'p'+tz,  4, False)
      Eof(False, sx*sz)
     
      # Horz floor of slot on top of cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  1, False)
      Prtv(tx+'m'+tz,  1, False)
      Prtv(tx+'m'+tz,  2, False)
      Prtv(tx+'p'+tz,  2, False)
      Eof(False, sx*sz)
      
      # Horz ceil of hole on main cube:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  5, False)
      Prtv(tx+'p'+tz,  6, False)
      Prtv(tx+'m'+tz,  6, False)
      Prtv(tx+'m'+tz,  5, False)
      Eof(False, sx*sz)
      
      # Top horz face of bar in hole:
      Bof(debug_faces)
      Prtv(tx+'p'+tz, 12, False)
      Prtv(tx+'p'+tz, 11, False)
      Prtv(tx+'m'+tz, 11, False)
      Prtv(tx+'m'+tz, 12, False)
      Eof(False, sx*sz)
      
      # Bot horz face of bar in hole:
      Bof(debug_faces)
      Prtv(tx+'p'+tz, 13, False)
      Prtv(tx+'p'+tz, 14, False)
      Prtv(tx+'m'+tz, 14, False)
      Prtv(tx+'m'+tz, 13, False)
      Eof(False, sx*sz)
        
      # Horz face on handle:
      Bof(debug_faces)
      Prtv(tx+'p'+tz,  7, False)
      Prtv(tx+'m'+tz,  7, False)
      Prtv(tx+'m'+tz,  8, False)
      Prtv(tx+'p'+tz,  8, False)
      Prtv(tx+'p'+tz,  7, False)
      Prtv(tx+'p'+tz,  9, True)
      Prtv(tx+'p'+tz, 10, False)
      Prtv(tx+'m'+tz, 10, False)
      Prtv(tx+'m'+tz,  9, False)
      Prtv(tx+'p'+tz,  9, False)
      Eof(True, sx*sz)
  
  sys.stderr.write("wrote %d faces (expected %d)\n" % (len(Fnrm)-1, Nf_exp))
  assert len(Fnrm) == Nf_exp + 1
  assert len(Fbar) == Nf_exp + 1

  sys.stderr.write("found %d non-ghost edges (expected %d)\n" % (len(E)-1, Ne_exp))
  assert len(E) == Ne_exp + 1

  sys.stderr.write("found %d ghost edges (expected %d)\n" % (len(G)-1, Ng_exp))
  assert len(G) == Ng_exp + 1

  return Fnrm, Fbar, E, G
  
def face_normal_and_barycenter(Fp):
  # Compute normal {nrm} and barycenter {bar}
  # given a list of vertices {Fp[0..deg-1]}
  deg = len(Fp)
  assert deg >= 3
  po = Fp[0][:3]
  pa = Fp[1][:3]
  pb = Fp[2][:3]
  u = rn.sub(pa,po)
  v = rn.sub(pb,po)
  nrm, sz = rn.dir(rn.cross3d(u, v))

  bar = (0, 0, 0)
  for jv in range(deg):
    pj = Fp[jv][:3]
    bar = rn.mix(1.0, bar, 1.0/deg, pj)
    
  # Check planarity:
  for jv in range(deg):
    pj = Fp[jv][:3]
    u = rn.sub(pj, bar)
    s = rn.dot(u, nrm)
    assert abs(s) < 0.5*eps, f"face is not planar s = {s}"
  return nrm, bar
  
def write_triangles_OBJ(wro, vtot, pock):
  # Writes the triangles of triangulations of the faces of the object to
  # OBJ file {wro}, assuming that the vertex indices are in the dict
  # {vtot}.
  
  assert not pock, "triangulaton with pockmarks not supported yet"
  
  # Returns a list {D} of all edges of all triangles, each 
  # as a pair {(kv_org,kv_dst)} with {kv_org < kv_dst}.
  # Element {D[0]} is not used.
  
  Nt_exp =  68 + 100 +  80   # Expected triangle count.
  Nd_exp = 3*Nt_exp // 2     # Expected triangulation edge count. 
  
  # Variables that are global to the nested functions:
  Ds = set()                # Triangulation edges.
  Tnrm = [ None ]           # Triangle normals.
  Tbar = [ None ]           # Triangle centers.
  debug_face = False 
  
  def Sved(p_org, p_dst):
    # Saves the edge from {p_org} to {p_dst} in {Ds},
    # if traversed in increasing vertex index sense.

    nonlocal Ds, Tnrm, Tbar, debug_face
    kv_org = p_org[3]
    kv_dst = p_dst[3]
    arc = (kv_org, kv_dst)
    if arc in Ds:
      sys.stderr.write(f"    !! duplicate arc {p_org[4]}--{p_dst[4]}\n")
      assert False
    Ds.add(arc)
    return None

  def Ptit(tit):
    # Title for a face or set of faces.

    nonlocal Ds, Tnrm, Tbar, debug_face
    wro.write("\n");
    wro.write("# %s\n" % tit);
    sys.stderr.write("writing %s \n" % (tit))
    return None
  
  def Prtt(key0, iv0, key1, iv1, key2, iv2, order):
    # Outputs triangle with vertices {vtot[key][iv]}
    # where {(key,iv)} is {(key0,iv0)}, {(key1,iv1)}, {(key2,iv2)}.
    # The vertices are listed in that order or the reverse
    # depending on whether {order} is {+1} or {-1}.

    nonlocal Ds, Tnrm, Tbar, debug_face
    # Save vertices in {Fp[0..2]} augmented with labels:
    Fp = []
    Fp.append(vtot[key0][iv0] + ( f"{key0}:{iv0:02d}", ))
    Fp.append(vtot[key1][iv1] + ( f"{key1}:{iv1:02d}", ))
    Fp.append(vtot[key2][iv2] + ( f"{key2}:{iv2:02d}", ))

    wro.write("f")
    for rv in range(3):
      jv = rv if order == +1 else 2-rv
      p = Fp[jv] 
      assert len(p) == 5
      kv = p[3]
      lab = p[4]
      if debug_face:
        sys.stderr.write(" v%04d = %s = ( %9.*f %9.*f %9.*f )\n" % (kv, lab, prec, p[0], prec, p[1], prec, p[2]));
      wro.write(" %d" % kv);
      Sved(Fp[jv], Fp[(jv + 3 + order) % 3])
    wro.write("\n");
    if debug_face: sys.stderr.write("\n")

    nrm, bar = face_normal_and_barycenter(Fp)
    Tnrm.append(nrm)
    Tbar.append(bar)
    return None
     
  # Here are the faces:
      
  Ptit(f"faces with normal ±X")
    
  for sx in +1, -1:
    tx = 'mop'[sx+1]

    Ptit(f" faces in halfspace sx = {sx:+d}")
    
    # Vert face of cube:
    Prtt(tx+'pm', 4,  tx+'pp', 4, tx+'pp', 7, sx)
    Prtt(tx+'mm', 4,  tx+'mm', 7, tx+'mp', 4, sx)
    Prtt(tx+'mm', 4,  tx+'pm', 4, tx+'pm', 7, sx)
    Prtt(tx+'mp', 7,  tx+'pp', 4, tx+'mp', 4, sx)
    Prtt(tx+'pp', 7,  tx+'pp', 4, tx+'mp', 7, sx)
    Prtt(tx+'pm', 4,  tx+'pp', 7, tx+'pm', 7, sx)
    Prtt(tx+'mm', 4,  tx+'pm', 7, tx+'mm', 7, sx)
    Prtt(tx+'mm', 7,  tx+'mp', 7, tx+'mp', 4, sx)
   
    # Centrad vert wall of hole on cube:
    Prtt(tx+'pp', 5,  tx+'mp', 5, tx+'mm', 5, sx)
    Prtt(tx+'pp', 5,  tx+'mm', 5, tx+'pm', 5, sx)
    
    # Distad vert wall of hole on cube:
    Prtt(tx+'mp', 6,  tx+'pp', 6, tx+'pm', 6, sx)
    Prtt(tx+'mp', 6,  tx+'pm', 6, tx+'mm', 6, sx)
    
    # Centrad vert wall of hole on handle:
    Prtt(tx+'pp', 9,  tx+'mp', 9, tx+'pm', 9, sx)
    Prtt(tx+'mm', 9,  tx+'pm', 9, tx+'mp', 9, sx)
    
    # Distad vert wall of hole on handle:
    Prtt(tx+'pp',10,  tx+'pm',10, tx+'mp',10, sx)
    Prtt(tx+'mp',10,  tx+'pm',10, tx+'mm',10, sx)
    
    # Extreme vert wall of handle:
    Prtt(tx+'pp', 8,  tx+'mp', 8, tx+'mm', 8, sx)
    Prtt(tx+'mm', 8,  tx+'pm', 8, tx+'pp', 8, sx)

    for sz in +1, -1:
      tz = 'mop'[sz+1]
    
      Ptit(f"    faces in quadrant sx = {sx:+d} sz = {sz:+d}")
      
      # Centrad vert wall of slot on face of cube:
      Prtt(tx+'m'+tz, 1,  tx+'p'+tz, 1,  tx+'m'+tz, 0, sx*sz)
      Prtt(tx+'p'+tz, 1,  tx+'p'+tz, 0,  tx+'m'+tz, 0, sx*sz)
      
      # Distad vert wall of slot on face of cube:
      debug_face = sx == -1 and sz == -1
      Prtt(tx+'p'+tz, 2,  tx+'m'+tz, 2,  tx+'m'+tz, 3, sx*sz)
      debug_face = False
      Prtt(tx+'p'+tz, 3,  tx+'p'+tz, 2,  tx+'m'+tz, 3, sx*sz)
      
      # Distad vert wall of bar in hole:
      Prtt(tx+'p'+tz,12,  tx+'m'+tz,12,  tx+'p'+tz,14, sx*sz)
      Prtt(tx+'p'+tz,14,  tx+'m'+tz,12,  tx+'m'+tz,14, sx*sz)
       
      # Centrad vert wall of bar in hole:
      Prtt(tx+'m'+tz,11,  tx+'p'+tz,11,  tx+'p'+tz,13, sx*sz)
      Prtt(tx+'m'+tz,13,  tx+'m'+tz,11,  tx+'p'+tz,13, sx*sz)
   
    
  Ptit(f"faces with normal ±Y")

  for sy in +1, -1:
    ty = 'mop'[sy+1]
    
    Ptit(f"  faces in halfspace sy = {sy:+d}")
    
    # Vert face of main cube:
    Ptit(f"    big complicated face")

    Prtt('m'+ty+'m', 4,  'm'+ty+'m', 2,  'm'+ty+'m', 3, sy)
    Prtt('m'+ty+'m', 4,  'm'+ty+'m', 6,  'm'+ty+'m', 2, sy)
    Prtt('m'+ty+'m', 4,  'm'+ty+'p', 6,  'm'+ty+'m', 6, sy)
    Prtt('m'+ty+'m', 4,  'm'+ty+'p', 4,  'm'+ty+'p', 6, sy)
    Prtt('m'+ty+'p', 6,  'm'+ty+'p', 4,  'm'+ty+'p', 2, sy)
    Prtt('m'+ty+'p', 2,  'm'+ty+'p', 4,  'm'+ty+'p', 3, sy)
                                      
    Prtt('m'+ty+'m', 2,  'm'+ty+'m', 5,  'm'+ty+'m', 1, sy)
    Prtt('m'+ty+'m', 2,  'm'+ty+'m', 6,  'm'+ty+'m', 5, sy)
    Prtt('m'+ty+'p', 6,  'm'+ty+'p', 2,  'm'+ty+'p', 5, sy)
    Prtt('m'+ty+'p', 5,  'm'+ty+'p', 2,  'm'+ty+'p', 1, sy)
                                      
    Prtt('m'+ty+'m', 0,  'm'+ty+'m', 1,  'o'+ty+'m', 0, sy)
    Prtt('o'+ty+'m', 0,  'm'+ty+'m', 1,  'o'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 1,  'm'+ty+'m', 1,  'm'+ty+'m', 5, sy)
    Prtt('o'+ty+'m', 1,  'm'+ty+'m', 5,  'm'+ty+'p', 5, sy)
    Prtt('o'+ty+'p', 0,  'm'+ty+'p', 5,  'm'+ty+'p', 1, sy)
    Prtt('o'+ty+'p', 0,  'm'+ty+'p', 1,  'm'+ty+'p', 0, sy)
                                      
    Prtt('m'+ty+'p', 5,  'o'+ty+'p', 0,  'p'+ty+'p', 5, sy)
    Prtt('p'+ty+'p', 5,  'o'+ty+'m', 1,  'm'+ty+'p', 5, sy)
                                      
    Prtt('p'+ty+'m', 0,  'o'+ty+'m', 0,  'p'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 0,  'o'+ty+'m', 1,  'p'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 1,  'p'+ty+'m', 5,  'p'+ty+'m', 1, sy)
    Prtt('o'+ty+'m', 1,  'p'+ty+'p', 5,  'p'+ty+'m', 5, sy)
    Prtt('o'+ty+'p', 0,  'p'+ty+'p', 1,  'p'+ty+'p', 5, sy)
    Prtt('o'+ty+'p', 0,  'p'+ty+'p', 0,  'p'+ty+'p', 1, sy)
                                      
    Prtt('p'+ty+'m', 1,  'p'+ty+'m', 5,  'p'+ty+'m', 2, sy)
    Prtt('p'+ty+'m', 2,  'p'+ty+'m', 5,  'p'+ty+'m', 6, sy)
    Prtt('p'+ty+'p', 5,  'p'+ty+'p', 2,  'p'+ty+'p', 6, sy)
    Prtt('p'+ty+'p', 2,  'p'+ty+'p', 5,  'p'+ty+'p', 1, sy)
                                      
    Prtt('p'+ty+'m', 4,  'p'+ty+'m', 3,  'p'+ty+'m', 2, sy)
    Prtt('p'+ty+'m', 4,  'p'+ty+'m', 2,  'p'+ty+'m', 6, sy)
    Prtt('p'+ty+'m', 4,  'p'+ty+'m', 6,  'p'+ty+'p', 6, sy)
    Prtt('p'+ty+'m', 4,  'p'+ty+'p', 6,  'p'+ty+'p', 4, sy)
    Prtt('p'+ty+'p', 4,  'p'+ty+'p', 6,  'p'+ty+'p', 2, sy)
    Prtt('p'+ty+'p', 4,  'p'+ty+'p', 2,  'p'+ty+'p', 3, sy)
    
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sx = {sx:+d} sy = {sy:+d}")
     
      # Vert outer face of handle:
      Prtt(tx+ty+'p', 7,  tx+ty+'p', 8,  tx+ty+'m', 8, sx*sy)
      Prtt(tx+ty+'p', 7,  tx+ty+'m', 8,  tx+ty+'m', 7, sx*sy)
                                      
      # Vert inner face of handle:
      Prtt(tx+ty+'p',10,  tx+ty+'p', 9,  tx+ty+'m', 9, sx*sy)
      Prtt(tx+ty+'p',10,  tx+ty+'m', 9,  tx+ty+'m',10, sx*sy)

      for sz in +1, -1:
        tz = 'mop'[sz+1]

        Ptit(f"      faces in octant sx = {sx:+d} sy = {sy:+d} sz = {sz:+d}")

        # Tip of bar in hole:
        Prtt(tx+ty+tz,13,  tx+ty+tz,11,  tx+ty+tz,14, sx*sy*sz)
        Prtt(tx+ty+tz,14,  tx+ty+tz,11,  tx+ty+tz,12, sx*sy*sz)
      
    
  Ptit(f"faces with normal ±Z")

  for sz in +1, -1:
    tz = 'mop'[sz+1]
    
    Ptit(f"  faces in halfspace sz = {sz:+d}")
    
    # Centrad horz face of main cube:
    Prtt('op'+tz, 0,  'om'+tz, 0,  'pp'+tz, 0, sz)
    Prtt('om'+tz, 0,  'pm'+tz, 0,  'pp'+tz, 0, sz)
    Prtt('om'+tz, 0,  'op'+tz, 0,  'mm'+tz, 0, sz)
    Prtt('op'+tz, 0,  'mp'+tz, 0,  'mm'+tz, 0, sz)
    
    for sx in +1, -1:
      tx = 'mop'[sx+1]
      
      Ptit(f"    faces in quadrant sx = {sx:+d} sz = {sz:+d}")

      # Distad horz face of main cube:
      Prtt(tx+'m'+tz, 3,  tx+'m'+tz, 4,  tx+'p'+tz, 4, sx*sz)
      Prtt(tx+'m'+tz, 3,  tx+'p'+tz, 4,  tx+'p'+tz, 3, sx*sz)
      
      # Horz floor face of slot on face of cube:
      Prtt(tx+'m'+tz, 1,  tx+'m'+tz, 2,  tx+'p'+tz, 2, sx*sz)
      Prtt(tx+'m'+tz, 1,  tx+'p'+tz, 2,  tx+'p'+tz, 1, sx*sz)
      
      # Distad floor of bar in hole:
      Prtt(tx+'p'+tz,12,  tx+'p'+tz,11,  tx+'m'+tz,12, sx*sz)
      Prtt(tx+'p'+tz,11,  tx+'m'+tz,11,  tx+'m'+tz,12, sx*sz)
       
      # Centrad floor of bar in hole:
      Prtt(tx+'p'+tz,13,  tx+'p'+tz,14,  tx+'m'+tz,13, sx*sz)
      Prtt(tx+'p'+tz,14,  tx+'m'+tz,14,  tx+'m'+tz,13, sx*sz)

      # Horz face of hole on main cube:
      Prtt(tx+'p'+tz, 5,  tx+'p'+tz, 6,  tx+'m'+tz, 5, sx*sz)
      Prtt(tx+'m'+tz, 5,  tx+'p'+tz, 6,  tx+'m'+tz, 6, sx*sz)
      
      # Horz face on handle:
      Prtt(tx+'p'+tz, 9,  tx+'p'+tz, 8,  tx+'p'+tz, 7, sx*sz)
      Prtt(tx+'p'+tz, 9,  tx+'p'+tz,10,  tx+'p'+tz, 8, sx*sz)
      Prtt(tx+'m'+tz, 8,  tx+'p'+tz, 8,  tx+'p'+tz,10, sx*sz)
      Prtt(tx+'m'+tz, 8,  tx+'p'+tz,10,  tx+'m'+tz,10, sx*sz)
      Prtt(tx+'m'+tz, 7,  tx+'m'+tz, 8,  tx+'m'+tz,10, sx*sz)
      Prtt(tx+'m'+tz, 7,  tx+'m'+tz,10,  tx+'m'+tz, 9, sx*sz)
      Prtt(tx+'m'+tz, 7,  tx+'m'+tz, 9,  tx+'p'+tz, 7, sx*sz)
      Prtt(tx+'m'+tz, 9,  tx+'p'+tz, 9,  tx+'p'+tz, 7, sx*sz)
    
  Nt_cmp = len(Tnrm) - 1
  sys.stderr.write("wrote %d triangles (expected %d)\n" % (Nt_cmp, Nt_exp))
  assert Nt_cmp == Nt_exp
    
  # Convert to list, leaving only one arc per edge:
  D = [ None ] + [ arc for arc in Ds if arc[0] < arc[1] ]
  
  Nd_cmp = len(D) - 1
  sys.stderr.write("wrote %d edges (expected %d)\n" % (Nd_cmp, Nd_exp))
  #assert Nd_cmp == Nd_exp

  return Tnrm, Tbar, D