code wiki / (root) / nx_gm_cloth_lib.nx
`nx_gm_cloth_lib.nx` source

↩ module page · 309 lines · 15224 B
1// nx_gm_cloth_lib.nx -- LIB: the PG29 cloth XPBD solver (gm_xpbd_step and its edge/area
2// rulers). Split out of nx_nxa_garment so the GATE can import it without colliding with the
3// generator's main(); the garment organ imports this lib and carries a REAL call (the
4// xpbdcheck verb), so its PG29 watch-row flip is backed by wiring, not a name in a comment.
5// license_tier: ORIGINAL   No hw writes (Rule 26).
6import "nx_softtissue.nx"
7// ---- PG29: THE CLOTH XPBD SOLVER (gm_xpbd_step) -- appended 2026-08-26 ----------------------
8// nx_nxa_garment GENERATED a fitted dress; NOTHING simulated it (the page's own comment:
9// "moves with the skin, zero physics", and GPIN is emitted but unconsumed). This is the solver
10// half of PG29, the organ symbol its watch row names. It is a LIB surface inside this organ:
11// main() above is untouched, so the generator's behaviour is bit-identical by construction.
12//
13// SHAPE: standard XPBD (Mueller et al., Position-Based Simulation of Compliant Constrained
14// Dynamics) -- predict under gravity, then iterate distance-constraint projections whose
15// stiffness enters ONLY through a compliance alpha, converted to alpha-tilde = alpha/dt^2 by
16// st_alpha_tilde -- COMPOSED from nx_softtissue, the estate's one XPBD lane, never re-derived.
17// All positions/velocities are i64 in the caller's mesh units; no float anywhere.
18//
19// NO MAGIC NUMBERS: gravity, dt, iteration count and both compliances are CALLER ARGUMENTS.
20// The solver owns only structure. Absolute fabric parameters (a cited bending stiffness) are a
21// NAMED follow-on -- gm_stiffness_mask, the second half of PG29 -- and until they are banked
22// the gate proves RELATIVE physics: a fold happens, area is conserved, stiffness changes the
23// fold. That is exactly the plan row's pre-declared done-rule.
24
25// ---- edge list of the GU x GV ring-grid, DERIVED from the same topology main() emits ----
26// classes: 0 = structural (ring + vertical), 1 = shear (quad diagonals), 2 = bending (2-hop
27// vertical). Rest lengths are measured from the rest verts handed in -- never assumed.
28// Row format per edge: [a, b, rest, class]. Returns the edge count written.
29const GM_EROW: i64 = 4
30const GM_C_STRUCT: i64 = 0
31const GM_C_SHEAR: i64 = 1
32const GM_C_BEND: i64 = 2
33// contact skin = the integer-truncation bound of the radial push (see gm_capsule_project)
34const GM_SKIN: i64 = 2
35
36func gm_isqrt(n: i64) -> i64 {
  if n <= 0 { return 0 }
  var x: i64 = n
  if x > 4096 { x = n/2 }
  var it: i64 = 0
  while it < 40 {
      if x > 0 { let y: i64 = (x + n/x)/2; if y > 0 { x = y } }
      it = it + 1
  }
  return x
46}
47
48func gm_dist(v: *i64, a: i64, b: i64) -> i64 {
  let dx: i64 = v[a*3] - v[b*3]
  let dy: i64 = v[a*3+1] - v[b*3+1]
  let dz: i64 = v[a*3+2] - v[b*3+2]
  return gm_isqrt(dx*dx + dy*dy + dz*dz)
53}
54
55// build the edge list over an nu x nv2 ring grid (column-major ring index u + row r, vertex
56// id = r*nu + u, ring wraps u -> (u+1)%nu). Caller supplies rest verts; rests are MEASURED.
57func gm_edges_build(verts: *i64, nu: i64, nv2: i64, edges: *i64) -> i64 {
  var ne: i64 = 0
  var r: i64 = 0
  while r < nv2 {
      var u: i64 = 0
      while u < nu {
          let a: i64 = r*nu + u
          let u1: i64 = (u + 1) % nu
          // ring structural
          let br: i64 = r*nu + u1
          edges[ne*GM_EROW] = a; edges[ne*GM_EROW+1] = br
          edges[ne*GM_EROW+2] = gm_dist(verts, a, br); edges[ne*GM_EROW+3] = GM_C_STRUCT
          ne = ne + 1
          if r + 1 < nv2 {
              // vertical structural
              let bv: i64 = (r+1)*nu + u
              edges[ne*GM_EROW] = a; edges[ne*GM_EROW+1] = bv
              edges[ne*GM_EROW+2] = gm_dist(verts, a, bv); edges[ne*GM_EROW+3] = GM_C_STRUCT
              ne = ne + 1
              // shear diagonals of the quad (a, ring-next, below, below-ring-next)
              let bd: i64 = (r+1)*nu + u1
              edges[ne*GM_EROW] = a; edges[ne*GM_EROW+1] = bd
              edges[ne*GM_EROW+2] = gm_dist(verts, a, bd); edges[ne*GM_EROW+3] = GM_C_SHEAR
              ne = ne + 1
              edges[ne*GM_EROW] = br; edges[ne*GM_EROW+1] = bv
              edges[ne*GM_EROW+2] = gm_dist(verts, br, bv); edges[ne*GM_EROW+3] = GM_C_SHEAR
              ne = ne + 1
          }
          if r + 2 < nv2 {
              // bending: 2-hop vertical resists the fold; its compliance is the stiffness knob
              let bb: i64 = (r+2)*nu + u
              edges[ne*GM_EROW] = a; edges[ne*GM_EROW+1] = bb
              edges[ne*GM_EROW+2] = gm_dist(verts, a, bb); edges[ne*GM_EROW+3] = GM_C_BEND
              ne = ne + 1
          }
          u = u + 1
      }
      r = r + 1
  }
  return ne
97}
98
99// ---- ONE XPBD STEP ---------------------------------------------------------------------------
100// verts/vel: i64 xyz triples, caller's mesh units. pins: vertex ids whose position is HELD
101// (infinite mass -- the GPIN contract). gy_per_tick2: gravity as velocity-delta per tick in mesh
102// units (sign carries direction; the garment's up-axis is +z, so a falling sheet passes negative).
103// alpha_q12_struct / alpha_q12_bend: XPBD compliance per class in the SAME Q12 the tissue lane
104// banks (st_alpha_tilde's own contract); shear uses the structural value. dt_us: microseconds per
105// tick, handed to st_alpha_tilde so stiffness is TIMESTEP-INDEPENDENT -- the whole point of XPBD.
106func gm_xpbd_step(verts: *i64, vel: *i64, nvv: i64,
                edges: *i64, ne: i64,
                pins: *i64, npins: i64,
                gz_per_tick: i64, dt_us: i64, iters: i64,
                alpha_q12_struct: i64, alpha_q12_bend: i64) -> i64 {
  // 1) predict: v += g; x += v (semi-implicit, the estate's integrator convention)
  var i: i64 = 0
  while i < nvv {
      vel[i*3+2] = vel[i*3+2] + gz_per_tick
      verts[i*3]   = verts[i*3]   + vel[i*3]
      verts[i*3+1] = verts[i*3+1] + vel[i*3+1]
      verts[i*3+2] = verts[i*3+2] + vel[i*3+2]
      i = i + 1
  }
  // bank pre-projection positions for the velocity update
  // (velocity = (x_new - x_pre_step)/1; x_pre_step = x_now - v after predict)
  // 2) pin restore BEFORE projection so constraints pull the cloth toward pins, not pins away
  var pp: i64 = 0
  while pp < npins {
      let pv: i64 = pins[pp]
      verts[pv*3]   = verts[pv*3]   - vel[pv*3]
      verts[pv*3+1] = verts[pv*3+1] - vel[pv*3+1]
      verts[pv*3+2] = verts[pv*3+2] - vel[pv*3+2]
      vel[pv*3] = 0; vel[pv*3+1] = 0; vel[pv*3+2] = 0
      pp = pp + 1
  }
  // 3) iterate distance projections. XPBD correction for equal unit masses with compliance:
  //    dl = (len - rest) / (2 + alpha_tilde)  applied +/- along the edge direction.
  //    alpha_tilde in Q12 via st_alpha_tilde (dsum = 2 unit masses), composed not re-derived.
  let at_s: i64 = st_alpha_tilde(alpha_q12_struct, 2, dt_us)
  let at_b: i64 = st_alpha_tilde(alpha_q12_bend, 2, dt_us)
  var it: i64 = 0
  while it < iters {
      var e: i64 = 0
      while e < ne {
          let a: i64 = edges[e*GM_EROW]
          let b: i64 = edges[e*GM_EROW+1]
          let rest: i64 = edges[e*GM_EROW+2]
          let cls: i64 = edges[e*GM_EROW+3]
          var atl: i64 = at_s
          if cls == GM_C_BEND { atl = at_b }
          let dx: i64 = verts[b*3] - verts[a*3]
          let dy: i64 = verts[b*3+1] - verts[a*3+1]
          let dz: i64 = verts[b*3+2] - verts[a*3+2]
          var ln: i64 = gm_isqrt(dx*dx + dy*dy + dz*dz)
          if ln < 1 { ln = 1 }
          // correction magnitude per endpoint, Q12 compliance in the denominator:
          // dl = (ln-rest)*4096 / ((2*4096) + atl) per endpoint pair -> split half-half below
          let dl2: i64 = (ln - rest)*4096/(8192 + atl)
          // apply half to each endpoint, along the edge
          let cxa: i64 = dx*dl2/ln
          let cya: i64 = dy*dl2/ln
          let cza: i64 = dz*dl2/ln
          var afree: i64 = 1
          var bfree: i64 = 1
          var q: i64 = 0
          while q < npins {
              if pins[q] == a { afree = 0 }
              if pins[q] == b { bfree = 0 }
              q = q + 1
          }
          if afree == 1 { if bfree == 1 {
              verts[a*3] = verts[a*3] + cxa/2; verts[a*3+1] = verts[a*3+1] + cya/2; verts[a*3+2] = verts[a*3+2] + cza/2
              verts[b*3] = verts[b*3] - cxa/2; verts[b*3+1] = verts[b*3+1] - cya/2; verts[b*3+2] = verts[b*3+2] - cza/2
          } }
          if afree == 1 { if bfree == 0 {
              verts[a*3] = verts[a*3] + cxa; verts[a*3+1] = verts[a*3+1] + cya; verts[a*3+2] = verts[a*3+2] + cza
          } }
          if afree == 0 { if bfree == 1 {
              verts[b*3] = verts[b*3] - cxa; verts[b*3+1] = verts[b*3+1] - cya; verts[b*3+2] = verts[b*3+2] - cza
          } }
          e = e + 1
      }
      it = it + 1
  }
  // 4) velocity update from actual displacement this tick: v = x - x_prev, where x_prev is
  // reconstructed as (x_after_predict - v_predict) per component -- we banked nothing, so
  // recompute: x_prev = x_now_would_need_banking. Instead XPBD's standard form: the predict
  // already wrote v; projections changed x without touching v, so fold the projection into v:
  // v += (x_projected - x_predicted). x_predicted is unavailable without a bank -- so this
  // solver banks nothing and instead DAMPS v toward the projected motion implicitly next tick
  // via the constraint pass. HONEST LIMIT, named: without the position bank the velocity keeps
  // the free-fall component; the gate's area/fold teeth measure POSITIONS, which are exact.
  // The bank costs nvv*3 words; gm_xpbd_step_banked below carries it and IS the recommended
  // entry point. This unbanked form exists so a memory-tight caller can still fold.
  return 0
192}
193
194// The recommended entry: identical, plus a caller-provided scratch bank (nvv*3 i64) so the
195// velocity update is the true XPBD one: v = x_final - x_start (per tick).
196func gm_xpbd_step_banked(verts: *i64, vel: *i64, nvv: i64,
                       edges: *i64, ne: i64,
                       pins: *i64, npins: i64,
                       gz_per_tick: i64, dt_us: i64, iters: i64,
                       alpha_q12_struct: i64, alpha_q12_bend: i64,
                       bank: *i64) -> i64 {
  var i: i64 = 0
  while i < nvv {
      bank[i*3] = verts[i*3]; bank[i*3+1] = verts[i*3+1]; bank[i*3+2] = verts[i*3+2]
      i = i + 1
  }
  gm_xpbd_step(verts, vel, nvv, edges, ne, pins, npins, gz_per_tick, dt_us, iters,
               alpha_q12_struct, alpha_q12_bend)
  i = 0
  while i < nvv {
      vel[i*3]   = verts[i*3]   - bank[i*3]
      vel[i*3+1] = verts[i*3+1] - bank[i*3+1]
      vel[i*3+2] = verts[i*3+2] - bank[i*3+2]
      i = i + 1
  }
  return 0
217}
218
219// total triangle area x2 of the ring grid (cross-product magnitudes), the conservation ruler
220// the gate's area tooth reads. Integer, exact up to isqrt rounding, and the SAME function
221// measures before and after so the comparison cannot drift between two rulers.
222func gm_area2(verts: *i64, nu: i64, nv2: i64) -> i64 {
  var total: i64 = 0
  var r: i64 = 0
  while r < nv2 - 1 {
      var u: i64 = 0
      while u < nu {
          let u1: i64 = (u + 1) % nu
          let a: i64 = r*nu + u
          let b: i64 = r*nu + u1
          let c: i64 = (r+1)*nu + u
          let d: i64 = (r+1)*nu + u1
          // quad = two triangles (a,b,c) and (b,d,c)
          var t9: i64 = 0
          while t9 < 2 {
              var p1: i64 = a; var p2: i64 = b; var p3: i64 = c
              if t9 == 1 { p1 = b; p2 = d; p3 = c }
              let e1x: i64 = verts[p2*3] - verts[p1*3]
              let e1y: i64 = verts[p2*3+1] - verts[p1*3+1]
              let e1z: i64 = verts[p2*3+2] - verts[p1*3+2]
              let e2x: i64 = verts[p3*3] - verts[p1*3]
              let e2y: i64 = verts[p3*3+1] - verts[p1*3+1]
              let e2z: i64 = verts[p3*3+2] - verts[p1*3+2]
              let cx: i64 = e1y*e2z - e1z*e2y
              let cy: i64 = e1z*e2x - e1x*e2z
              let cz: i64 = e1x*e2y - e1y*e2x
              total = total + gm_isqrt(cx*cx + cy*cy + cz*cz)
              t9 = t9 + 1
          }
          u = u + 1
      }
      r = r + 1
  }
  return total
255}
256
257// ---- CLOTH-OVER-BODY CONTACT: capsule projection --------------------------------------------
258// The world page's body is a skinned figure whose limbs and torso a CAPSULE approximates -- the
259// same primitive st_set_touch's sphere probe already uses one rung down. Projection, not force:
260// after the constraint iterations, any cloth vert inside the capsule is pushed to its surface
261// along the radial, exactly the PBD contact convention. Called once per tick per capsule by the
262// stepper's caller (capsule count is the CALLER's anatomy, not this lib's constant).
263// Capsule: segment (ax,ay,az)-(bx,by,bz) with radius rad, all in the caller's mesh units.
264// Returns the number of verts projected -- 0 on a clean pass, so a gate can assert BOTH ways.
265func gm_capsule_project(verts: *i64, vel: *i64, nvv: i64,
                      ax: i64, ay: i64, az: i64, bx: i64, by: i64, bz: i64, rad: i64) -> i64 {
  let ux: i64 = bx - ax
  let uy: i64 = by - ay
  let uz: i64 = bz - az
  var uu: i64 = ux*ux + uy*uy + uz*uz
  if uu < 1 { uu = 1 }
  var hits: i64 = 0
  var i: i64 = 0
  while i < nvv {
      let px: i64 = verts[i*3] - ax
      let py: i64 = verts[i*3+1] - ay
      let pz: i64 = verts[i*3+2] - az
      // closest point on the segment: t in [0,1] scaled by 4096 (Q12, matching the lib's
      // compliance quantum so the lib carries ONE fixed-point convention)
      var tq: i64 = (px*ux + py*uy + pz*uz)*4096/uu
      if tq < 0 { tq = 0 }
      if tq > 4096 { tq = 4096 }
      let cx: i64 = ax + ux*tq/4096
      let cy: i64 = ay + uy*tq/4096
      let cz: i64 = az + uz*tq/4096
      let dx: i64 = verts[i*3] - cx
      let dy: i64 = verts[i*3+1] - cy
      let dz: i64 = verts[i*3+2] - cz
      let d2: i64 = dx*dx + dy*dy + dz*dz
      if d2 < rad*rad {
          var d: i64 = gm_isqrt(d2)
          if d < 1 { d = 1 }
          // push to the surface along the radial; kill the inward velocity component the same
          // way the XPBD clamp does (a projected vert must not re-enter next tick on inertia).
          // GM_SKIN: integer division truncates toward zero, so each axis loses under 1 unit and
          // the radial lands up to ~2 units INSIDE the surface -- measured as exactly the T7
          // failure. The skin is that truncation bound, DERIVED not tuned: push to rad+GM_SKIN
          // so the truncated result still rests at-or-outside rad.
          let rs: i64 = rad + GM_SKIN
          verts[i*3]   = cx + dx*rs/d
          verts[i*3+1] = cy + dy*rs/d
          verts[i*3+2] = cz + dz*rs/d
          vel[i*3] = 0; vel[i*3+1] = 0; vel[i*3+2] = 0
          hits = hits + 1
      }
      i = i + 1
  }
  return hits
309}
nx_gm_cloth_lib.nx source

`nx_gm_cloth_lib.nx` source
