diff --git a/aephysics/manifold/module.ae b/aephysics/manifold/module.ae index b8c5eff..cd55a6a 100644 --- a/aephysics/manifold/module.ae +++ b/aephysics/manifold/module.ae @@ -222,14 +222,15 @@ clip_points() -> ptr { return native.scratch(g_clip_points, MAX_CLIP_POINTS * si // A's edges, gathered once before the edge-pair loop reads them a // hundred times (see edge_axes). var g_edge_axes: ptr = null -edge_axes() -> ptr { return native.scratch(g_edge_axes, AXIS_FIELDS * AXIS_SLOTS * 4) } +edge_axes() -> ptr { return native.scratch(g_edge_axes, EDGE_SCRATCH_BYTES) } +b_frame_scratch() -> ptr { return edge_axes() + AXIS_FIELDS * AXIS_SLOTS * 4 } scratch_ready() { if g_clip1 == null { g_clip1 = basics.alloc(native.MAX_WORKERS * (MAX_CLIP_POINTS * sizeof(ClipVertex))) g_clip2 = basics.alloc(native.MAX_WORKERS * (MAX_CLIP_POINTS * sizeof(ClipVertex))) g_clip_points = basics.alloc(native.MAX_WORKERS * (MAX_CLIP_POINTS * sizeof(LocalManifoldPoint))) - g_edge_axes = basics.alloc(native.MAX_WORKERS * (AXIS_FIELDS * AXIS_SLOTS * 4)) + g_edge_axes = basics.alloc(native.MAX_WORKERS * EDGE_SCRATCH_BYTES) } } @@ -255,24 +256,32 @@ reserve_scratch() { scratch_ready() } // edges at a time (std.lanes, Aether 0.706). Each lane computes what the // edge computed on its own, in the same order, so the axis a pair // chooses is the axis it chose before: the tests hold the traces. -const AXIS_SLOTS = 64 // MAX_HULL_EDGES / 2, already a multiple of four +// A hull has up to MAX_HULL_EDGES full edges (2 * MAX_HULL_EDGES half-edges, +// the builder's limit), so the table holds that many, a multiple of four. +const AXIS_SLOTS = 128 const AXIS_FIELDS = 14 // Where each array starts, in elements: x, y, z of the direction, of // both face normals and of the origin, then the tolerance. -const AXIS_DIR_X = 0 -const AXIS_DIR_Y = 64 -const AXIS_DIR_Z = 128 -const AXIS_N0_X = 192 -const AXIS_N0_Y = 256 -const AXIS_N0_Z = 320 -const AXIS_N1_X = 384 -const AXIS_N1_Y = 448 -const AXIS_N1_Z = 512 -const AXIS_V0_X = 576 -const AXIS_V0_Y = 640 -const AXIS_V0_Z = 704 -const AXIS_TOL = 768 -const AXIS_EDGE = 832 // the candidate's half-edge index (exact in an f32) +const AXIS_DIR_X = 0 * AXIS_SLOTS +const AXIS_DIR_Y = 1 * AXIS_SLOTS +const AXIS_DIR_Z = 2 * AXIS_SLOTS +const AXIS_N0_X = 3 * AXIS_SLOTS +const AXIS_N0_Y = 4 * AXIS_SLOTS +const AXIS_N0_Z = 5 * AXIS_SLOTS +const AXIS_N1_X = 6 * AXIS_SLOTS +const AXIS_N1_Y = 7 * AXIS_SLOTS +const AXIS_N1_Z = 8 * AXIS_SLOTS +const AXIS_V0_X = 9 * AXIS_SLOTS +const AXIS_V0_Y = 10 * AXIS_SLOTS +const AXIS_V0_Z = 11 * AXIS_SLOTS +const AXIS_TOL = 12 * AXIS_SLOTS +const AXIS_EDGE = 13 * AXIS_SLOTS // the candidate's half-edge index (exact in an f32) +// The table, then room for B's face normals and vertices taken into A's frame +// after it: hull.MAX_HULL_FACES + hull.MAX_HULL_VERTICES of them, written out +// (a constant naming another module's is expanded where this module is +// imported, and not every importer imports the hull). +const B_FRAME_SLOTS = 256 +const EDGE_SCRATCH_BYTES = AXIS_FIELDS * AXIS_SLOTS * 4 + B_FRAME_SLOTS * 24 // --- the incident face and clipping ---------------------------------------------------- @@ -1093,11 +1102,10 @@ compute_separating_axis(hull_a: *HullData, hull_b: *HullData, xf_b: Transform, e f32 eps32 = eps if axis_count == 0 { return res } // B's face normals and vertices into A's frame, negated, once (Box3D - // adec25b) rather than four products for every B edge; they live in - // the first clip polygon's scratch, idle until the contact is built - // (a hull's faces and vertices number its edges plus two, at most 66). + // adec25b) rather than four products for every B edge, into the room + // after the edge table. scratch_ready() - b_in_a = clip1() as Vec3[] + b_in_a = b_frame_scratch() as Vec3[] k = 0 while k < hull_b.face_count { b_in_a[k] = math.neg(math.mul_mv(r, planes_b[k].normal)) diff --git a/tests/test_manifold.ae b/tests/test_manifold.ae index d62b548..102355f 100644 --- a/tests/test_manifold.ae +++ b/tests/test_manifold.ae @@ -676,6 +676,88 @@ test_edge_axis_oracle() { hull.destroy_hull(a) } +// A 32-sided cylinder hull: 64 vertices, 34 faces, 96 edges, past the 64 +// the edge-axis table used to hold (its later edges wrote over the earlier +// ones' fields and past the worker's block). +cylinder_hull(radius: float, half_height: float) -> *HullData { + block = calloc(64, sizeof(Vec3)) + p = block as Vec3[] + i = 0 + while i < 32 { + angle = 2.0 * math.PI * (i as float) / 32.0 + x = radius * cos(angle) + z = radius * sin(angle) + p[2 * i] = math.vec3(x, 0.0 - half_height, z) + p[2 * i + 1] = math.vec3(x, half_height, z) + i = i + 1 + } + h = hull.create_hull(p, 64, 64) + free(block) + return h +} + +// The best edge-pair separation by brute force, with the SAT's own rule +// for an axis (arcs crossing past its epsilon, not near parallel). +brute_edge_separation(a: *HullData, b: *HullData, xf: Transform) -> float { + r = math.make_matrix_from_quat(xf.q) + planes_a = hull.hull_planes(a) + points_a = hull.hull_points(a) + edges_a = hull.hull_edges(a) + planes_b = hull.hull_planes(b) + points_b = hull.hull_points(b) + edges_b = hull.hull_edges(b) + tol = manifold.PARALLEL_EDGE_TOL * manifold.PARALLEL_EDGE_TOL + eps = 0.0 - 0.0001 + best = 0.0 - math.MAX_FLOAT + i = 0 + while i < a.edge_count { + av0 = points_a[edges_a[i].origin] + ad = math.sub(points_a[edges_a[i + 1].origin], av0) + n0 = planes_a[edges_a[i].face].normal + n1 = planes_a[edges_a[i + 1].face].normal + j = 0 + while j < b.edge_count { + c = math.neg(math.mul_mv(r, planes_b[edges_b[j].face].normal)) + d = math.neg(math.mul_mv(r, planes_b[edges_b[j + 1].face].normal)) + bv0 = math.neg(math.add(math.mul_mv(r, points_b[edges_b[j].origin]), xf.p)) + bv1 = math.neg(math.add(math.mul_mv(r, points_b[edges_b[j + 1].origin]), xf.p)) + dc = math.sub(bv1, bv0) + cba = math.dot(c, ad) + dba = math.dot(d, ad) + adc = math.dot(n0, dc) + bdc = math.dot(n1, dc) + if cba * dba < eps && adc * bdc < eps && cba * bdc < eps && math.max_float(cba * cba, dba * dba) > tol * math.length_squared(ad) { + t = (0.0 - cba) / (dba - cba) + n = math.normalize(math.add(c, math.mul_sv(t, math.sub(d, c)))) + best = math.max_float(best, 0.0 - math.dot(math.add(av0, bv0), n)) + } + j = j + 2 + } + i = i + 2 + } + return best +} + +test_many_edges() { + a = cylinder_hull(0.5, 0.4) + b = cylinder_hull(0.4, 0.6) + ensure("many edges: hulls of 96 edges (${a.edge_count / 2}, ${b.edge_count / 2})", a != null && b != null && a.edge_count == 192 && b.edge_count == 192) + seed = 97531 + agree = 0 + trials = 300 + i = 0 + while i < trials { + xf = Transform { p: math.mul_sv(next_float(0.6, 1.4), next_direction()), q: exact_quat(next_direction(), next_float(0.1, 3.0)) } + query = manifold.compute_separating_axis(a, b, xf, false) + brute = brute_edge_separation(a, b, xf) + if math.abs_float(query.edge.separation - brute) < 0.0001 { agree = agree + 1 } + i = i + 1 + } + ensure("many edges: the SAT's edge phase agrees with brute force (${agree} of ${trials})", agree >= trials - 3) + hull.destroy_hull(a) + hull.destroy_hull(b) +} + // --- hull against capsule and sphere ----------------------------------------------------- // A thin capsule stabbed through the +x +y edge of a box: the isolated @@ -885,6 +967,7 @@ main() { test_overlap_never_empty() test_ridge_crossing() test_edge_axis_oracle() + test_many_edges() test_hull_capsule_edge_deep() test_sphere_hull_seam() test_rounded_pairs()