Merge branch 'main' into twidmer/more-contact-improvements

Signed-off-by: nvtw <110816143+nvtw@users.noreply.github.com>

Author: nvtw

CHANGELOG.md

@@ -17,6 +17,7 @@
 - Add `ViewerUSD(points_as_spheres=...)` to render `log_points` particles as a `UsdGeom.PointInstancer` of sphere prototypes; enabled by default (opt out with `points_as_spheres=False` for flat `UsdGeom.Points` splats)
 - Add experimental `newton.controllers` module with `ControllerBase` base class, `ControllerJointImpedance`, and `ControllerJointImpedanceModelFree` for GPU-accelerated, vectorized joint-space impedance control.
 - Add list-of-pattern and explicit-index selectors to `ArticulationView`.
+- Add full-surface (edge/face) rigid-soft contacts to `SolverVBD` proxy-body coupling under `SolverCoupledProxy` through shared or proxy-local collision pipelines. Proxy-particle coupling rejects full-surface contacts.
 - Add `newton[onnx]` for ONNX policy inference through Warp-NN; `ControllerNeuralMLP`, `ControllerNeuralLSTM`, and RL policy examples can run exported `.onnx` policies without requiring PyTorch for ONNX execution.
 - Add three VBD contact examples — `vbd_rigid_rigid_contact`, `vbd_soft_rigid_contact`, and `vbd_soft_rigid_mix_contact` — demonstrating rigid-rigid, soft (particle-rigid), and mixed cloth-bag contacts
 - Add masked rigid-body reset support to `SolverVBD`. (#3256)
@@ -47,6 +48,9 @@
 
 ### Changed
 
+- Decide collider visibility from USD `purpose` and visibility rather than from a bound render material. A collider whose `purpose` resolves to `default` is viewport geometry and is drawn; mark it `guide` to state that it is collision-only. Previously an unrelated visual elsewhere in the scene could make a collider vanish. `force_show_colliders` and `hide_collision_shapes` are unchanged.
+- Filter shared full-surface soft contacts per `SolverCoupled` entry, preserving them for capable solvers and dropping them for particle-only solvers or records spanning entries.
+- Require `warp-lang>=1.16.0`; upgrade Warp to version 1.16.0 or later.
 - Follow USD viewport `purpose` and visibility for imported visual geometry and colliders rather than inferring visibility from a bound render material. Visual shapes and Gaussian splats are drawn only for `default` and `proxy`; a collider whose `purpose` resolves to `default` is drawn, while `guide` identifies collision-only geometry. `force_show_colliders` and `hide_collision_shapes` are unchanged. (#3404, #3712)
 - Require `warp-lang>=1.16.0`; upgrade Warp to version 1.16.0 or later. (#3780)
 - Disable the implicit positive Dahl-friction defaults in `SolverVBD.register_custom_attributes()` (deprecated in 1.3.0): `vbd:dahl_eps_max` and `vbd:dahl_tau` now default to zero, and Dahl cable friction is enabled only where both are authored positive. Pass `dahl_defaults_enabled=True` to temporarily restore the old defaults; the compatibility mode will be removed in a future release.
@@ -57,7 +61,7 @@
 - Upgrade `mujoco` and `mujoco-warp` to 3.11.0. (#3725)
 - Optimize raycast/raytrace queries by restructuring ray-shape intersection into local-space primitives and compile specialized depth/shadow variants that skip unused surface-normal work (mesh shadows also use any-hit queries).
 - Speed up `SolverKamino` DVI contact coloring for dense manifolds requiring more than 64 colors.
-- Change experimental `SolverVBD` cable constraint slots from `[STRETCH=0, BEND=1]` to `[STRETCH=0, SHEAR=1, BEND=2, TWIST=3]`, allowing each stiffness and constraint mode to be configured independently. Existing cable calls using raw `slot=1` or `JointSlot.ANGULAR` now select shear; use `JointSlot.BEND` (now slot 2) to select bending.
+- Change experimental `SolverVBD` cable constraint slots from `[STRETCH=0, BEND=1]` to `[STRETCH=0, SHEAR=1, BEND=2, TWIST=3]`, allowing each stiffness and constraint mode to be configured independently while preserving the pre-split world-form parent solve Hessian for the common isotropic component of stretch/shear elasticity. Existing cable calls using raw `slot=1` or `JointSlot.ANGULAR` now select shear; use `JointSlot.BEND` (now slot 2) to select bending.
 - Map `shape_material_kf` to per-contact MuJoCo `solreffriction` in `SolverMuJoCo` (elliptic friction cones with Newton contacts); resolve `kf` with priority/`solmix`, treat a resolved `kf = 0` as frictionless, and use native MuJoCo contacts or a pyramidal cone to preserve the previous solref-inherited friction.
 - Load visual-only USD geometry outside rigid-body hierarchies as static shapes by default; pass `load_static_visual_shapes=False` to retain the previous body-associated-visuals-only behavior.
 - Speed up `Mesh.create_heightfield()` and `Mesh.create_terrain()` by building the vertex and index buffers in place, substantially reducing construction time and peak memory for large terrain grids such as those used by Isaac Lab.
@@ -95,6 +99,7 @@
 - Make deterministic collision pipelines cover hydroelastic contact generation and reduction, including unique reduced-contact sort keys and overflow-safe fixed-point pressure accumulation. (#3661)
 - Fix `SolverMuJoCo` retaining an invalid external-contact cache when its first step is captured in a CUDA graph. (#3768; fixes #3767)
 - Preserve box-box face contact manifolds under sub-microradian solver drift. (#3776)
+- Fix `SolverMuJoCo` overflowing MuJoCo's signed 32-bit collision masks when graph coloring requires the highest supported color; all 32 mask bits are now used before falling back to default collision masks for additional colors.
 - Convert `newton:mimicCoef0` from degrees to radians when the mimic follower joint is angular. Assets authored against the old behavior need the value rescaled to degrees.
 - Complete Kamino RCM traversal for large and disconnected systems and reuse the resulting permutation by default; set `reuse_permutation=False` to recompute it for changing matrix topology.
 - Bound Kamino DVI contact allocation with a per-world geometry heuristic instead of sizing every contact pair simultaneously; set `collision_detector.max_contacts_per_world` to override the inferred capacity.
@@ -106,6 +111,7 @@
 - Fix panel-parallel RCM-blocked LLT factorization hanging when a matrix ends in a partial tile.
 - Fix USD capsule, cylinder, and cone visual and site scaling to follow the authored primitive axis.
 - Fix MJCF contact pairs ignoring properties inherited from pair default classes.
+- Fix disabled USD colliders participating in particle collisions when visual shape loading is disabled.
 - Fix `ArticulationView.is_fixed_base` for roots with zero effective degrees of freedom, including fully locked D6 joints. (#3727)
 - Fix USD plane visual width and length to scale along the axes defined by the `UsdGeomPlane` schema, and orient X- and Y-axis plane visuals along the authored axis.
 - Validate `ArticulationView` mask shapes and devices before launching selection kernels. (#3448)

docs/concepts/coupling.rst

@@ -379,6 +379,11 @@ The coupled-solver framework is useful today, but it is still experimental:
   custom hook is available.
 - USD ownership, automatic coupled-solver construction, and high-level tuning
   guidance are not part of the experimental public API yet.
+- Full-surface (edge/face) rigid-soft contacts are consumed by
+  :class:`~newton.solvers.SolverVBD` only. The per-entry contact filter drops
+  them for other sub-solvers while preserving particle contacts, and keeps them
+  for VBD only when the entry owns every referenced corner. A record spanning
+  two entries is dropped by both. Proxy particles do not support them at all.
 
 Treat coupled solvers as an advanced feature for controlled experiments and
 solver integration work. Prefer focused regression tests and explicit scene

newton/_src/solvers/coupled/interface.py

@@ -44,6 +44,9 @@ def coupling_notify_input_state_update(state, flags, *, iteration_restart=False,
     def coupling_supports_inertial_property_refresh() -> bool: ...
 
 
+    def coupling_supports_full_surface_soft_contacts() -> bool: ...
+
+
     def coupling_rewind_proxy_body(
         body_local_to_proxy_global, state, coupling_forces, body_gravity_acceleration, dt
     ) -> None: ...
@@ -252,6 +255,10 @@ def coupling_supports_inertial_property_refresh(self) -> bool:
         """
         return False
 
+    def coupling_supports_full_surface_soft_contacts(self) -> bool:
+        """Return whether the solver consumes edge and face soft contacts."""
+        return False
+
     def coupling_eval_gravity_acceleration(
         self,
         out_body_acceleration: wp.array[wp.vec3] | None,

newton/_src/solvers/coupled/solver_coupled.py

@@ -2711,6 +2711,7 @@ def _contacts_for_entry(self, entry: SolverEntry, contacts: Contacts | None) ->
             return contacts
 
         filtered = self._ensure_entry_contact_buffer(entry, contacts)
+        keep_full_surface_contacts = filtered._enable_rigid_soft_full_surface_contact
         force_contact_update = int(self._entry_contact_sources.get(entry.name) is not contacts)
         if force_contact_update:
             self._entry_contact_sources[entry.name] = contacts
@@ -2849,6 +2850,7 @@ def _contacts_for_entry(self, entry: SolverEntry, contacts: Contacts | None) ->
                     entry.view.particle_flags,
                     int(ShapeFlags.COLLIDE_PARTICLES),
                     int(ParticleFlags.ACTIVE),
+                    int(keep_full_surface_contacts),
                     filtered.soft_contact_count,
                     filtered.soft_contact_particle,
                     filtered.soft_contact_shape,
@@ -2917,6 +2919,10 @@ def _ensure_entry_contact_buffer(self, entry: SolverEntry, contacts: Contacts) -
                 device=contacts.device,
             )
         filtered._contact_matching_mode = contacts.contact_matching_mode
+        filtered._enable_rigid_soft_full_surface_contact = bool(
+            contacts._enable_rigid_soft_full_surface_contact
+            and entry.solver.coupling_supports_full_surface_soft_contacts()
+        )
         return filtered
 
     @staticmethod
@@ -3674,6 +3680,7 @@ def _filter_soft_contacts_global_shape_ids_kernel(
     particle_flags: wp.array[wp.int32],
     collide_particles_mask: int,
     active_particle_mask: int,
+    keep_full_surface_contacts: int,
     dst_count: wp.array[wp.int32],
     dst_particle: wp.array[int],
     dst_shape: wp.array[int],
@@ -3694,14 +3701,24 @@ def _filter_soft_contacts_global_shape_ids_kernel(
 
     particle = src_particle[contact_id]
     shape = src_shape[contact_id]
-    if particle < 0 or shape < 0:
+    if particle < 0 and keep_full_surface_contacts == 0:
         return
-    if particle >= particle_flags.shape[0] or shape >= shape_flags.shape[0]:
+    if shape < 0 or shape >= shape_flags.shape[0]:
         return
     if (shape_flags[shape] & collide_particles_mask) == 0:
         return
-    if (particle_flags[particle] & active_particle_mask) == 0:
+
+    # Keep records only when the entry owns every referenced corner.
+    corners = src_indices[contact_id]
+    if corners[0] < 0:
         return
+    for i in range(3):
+        corner = corners[i]
+        if corner >= 0:
+            if corner >= particle_flags.shape[0]:
+                return
+            if (particle_flags[corner] & active_particle_mask) == 0:
+                return
 
     dst_id = wp.atomic_add(dst_count, 0, wp.int32(1))
     src_to_dst[contact_id] = dst_id

newton/_src/solvers/mujoco/solver_mujoco.py

@@ -5126,10 +5126,8 @@ def _color_collision_shapes(
                 graph_edge_indices=wp.array(graph_edges, dtype=wp.int32),
                 balance_colors=False,
             )
-            num_colors = 0
-            for group in color_groups:
-                num_colors += 1
-                shape_color[selected_shapes[group]] = num_colors
+            for color, group in enumerate(color_groups):
+                shape_color[selected_shapes[group]] = color
             if visualize_graph:
                 plot_graph(
                     vertices=np.arange(num_shapes),
@@ -5140,6 +5138,21 @@ def _color_collision_shapes(
 
         return shape_color
 
+    @staticmethod
+    def _collision_color_masks(color: int) -> tuple[int, int]:
+        """Encode a collision color as signed MuJoCo ``contype`` and ``conaffinity`` masks."""
+        if color >= 32:
+            return 1, 1
+
+        uint32_mask = (1 << 32) - 1
+        contype = 1 << color
+        conaffinity = uint32_mask & ~contype
+
+        def as_signed_int32(mask: int) -> int:
+            return mask if mask < (1 << 31) else mask - (1 << 32)
+
+        return as_signed_int32(contype), as_signed_int32(conaffinity)
+
     def get_max_contact_count(self) -> int:
         """Return the maximum number of rigid contacts that can be generated by MuJoCo."""
         if self.use_mujoco_cpu:
@@ -5613,9 +5626,6 @@ def joint_has_raw_limit_solref(dof_idx: int) -> bool:
         mimic_enabled = model.constraint_mimic_enabled.numpy()
         mimic_world = model.constraint_mimic_world.numpy()
 
-        INT32_MAX = np.iinfo(np.int32).max
-        collision_mask_everything = INT32_MAX
-
         # mapping from joint axis to actuator index
         # axis_to_actuator[i, 0] = position actuator index
         # axis_to_actuator[i, 1] = velocity actuator index
@@ -6131,11 +6141,9 @@ def add_geoms(newton_body_id: int):
                     geom_params["conaffinity"] = mujoco_mask_to_signed(int(compiled_collision_affinity[shape]))
                 else:
                     color = shape_color[shape]
-                    if color < 32:
-                        contype = 1 << color
-                        geom_params["contype"] = contype
-                        # collide with anything except shapes from the same color
-                        geom_params["conaffinity"] = collision_mask_everything & ~contype
+                    contype, conaffinity = self._collision_color_masks(int(color))
+                    geom_params["contype"] = contype
+                    geom_params["conaffinity"] = conaffinity
 
                 # set friction from Newton shape materials
                 mu = shape_mu[shape]

newton/_src/solvers/vbd/rigid_vbd_kernels.py

@@ -1088,8 +1088,10 @@ def evaluate_cable_stretch_shear_force_hessian(
     u = wp.quat_rotate_inv(q_wp, C_vec)
     psi = wp.cw_mul(k_diag, u) - C0_force_local + lambda_local
 
-    h_s = k_diag[0]
-    h_z = k_diag[2]
+    k_s = k_diag[0]
+    k_z = k_diag[2]
+    h_s = k_s
+    h_z = k_z
     if damping_active:
         inv_dt = 1.0 / dt
         x_p_prev = wp.transform_get_translation(X_wp_prev)
@@ -1103,11 +1105,24 @@ def evaluate_cable_stretch_shear_force_hessian(
     force = f_world if is_parent else -f_world
 
     t = _quat_rotate_local_z(q_wp)
-    K_eff = h_s * wp.identity(3, float) + (h_z - h_s) * wp.outer(t, t)
-    rx = wp.skew(r)
+    identity = wp.identity(3, float)
+    K_eff = h_s * identity + (h_z - h_s) * wp.outer(t, t)
     H_ll = K_eff
-    H_al = rx * K_eff
-    H_aa = wp.transpose(rx) * K_eff * rx
+    if is_parent:
+        # Isotropic elastic energy is frame-invariant, so it takes the parent-anchor arm that
+        # evaluate_linear_constraint_force_hessian also uses; min() extracts the largest such block
+        # leaving a PSD remainder. Damping keeps the material arm: u_prev is frozen one step back.
+        k_iso = wp.min(k_s, k_z)
+        K_material = K_eff - k_iso * identity
+        r_elastic = x_p - com_w
+        rx_material = wp.skew(r)
+        H_al = k_iso * wp.skew(r_elastic) + rx_material * K_material
+        H_aa = k_iso * (wp.length_sq(r_elastic) * identity - wp.outer(r_elastic, r_elastic))
+        H_aa = H_aa + wp.transpose(rx_material) * K_material * rx_material
+    else:
+        rx = wp.skew(r)
+        H_al = rx * K_eff
+        H_aa = wp.transpose(rx) * K_eff * rx
 
     torque = wp.cross(r, force)
     return force, torque, H_ll, H_al, H_aa

newton/_src/solvers/vbd/solver_vbd.py

@@ -868,6 +868,10 @@ def notify_model_changed(self, flags: ModelFlags | int) -> None:
     def coupling_supports_inertial_property_refresh(self) -> bool:
         return True
 
+    @override
+    def coupling_supports_full_surface_soft_contacts(self) -> bool:
+        return True
+
     def coupling_notify_input_state_update(
         self,
         state: State,
@@ -926,19 +930,6 @@ def coupling_prepare_proxy_contacts(
         contacts_freshly_detected: bool = False,
     ) -> Contacts | None:
         """Update rigid history cadence for proxy contacts."""
-        # Full-surface (edge/face) rigid-soft contacts are not yet harvested onto proxy particles:
-        # the proxy contact-force kernels consume only per-particle records (particle >= 0), so a soft
-        # edge/face contact's reaction on a proxy-coupled rigid body would be silently dropped. Fail
-        # loud until proxy harvesting consumes the unified records. Standalone SolverVBD (no proxy
-        # coupling) never reaches this hook and does support full-surface via the per-body path.
-        if contacts is not None and getattr(contacts, "_enable_rigid_soft_full_surface_contact", False):
-            raise NotImplementedError(
-                "Full-surface (edge/face) rigid-soft contacts are not yet supported with VBD proxy-particle "
-                "coupling (SolverCoupledProxy): the proxy contact-force harvest consumes only per-particle "
-                "records, so edge/face force feedback to proxy-coupled rigid bodies would be silently dropped. "
-                "Set enable_rigid_soft_full_surface_contact=False for the coupled proxy solve, or drive the "
-                "rigid bodies with standalone SolverVBD (which supports full-surface via the per-body path)."
-            )
         # Do not call super(); we can keep proxy-proxy collisions as we
         # are using a custom force harvesting hook
         self.set_rigid_history_update(bool(contacts_freshly_detected))
@@ -1032,7 +1023,8 @@ def coupling_harvest_proxy_wrenches(
                     self.body_particle_contact_material_kd,
                     self.body_particle_contact_material_mu,
                     contacts.soft_contact_count,
-                    contacts.soft_contact_particle,
+                    contacts.soft_contact_indices,
+                    contacts.soft_contact_barycentric,
                     contacts.soft_contact_shape,
                     contacts.soft_contact_body_pos,
                     contacts.soft_contact_body_vel,
@@ -1068,6 +1060,9 @@ def coupling_harvest_proxy_particle_forces(
         if self.model.particle_count == 0 or particle_local_to_proxy_global.shape[0] == 0:
             return
 
+        if contacts is not None:
+            contacts._assert_particle_only_soft_contacts("SolverVBD proxy-particle coupling")
+
         if (
             contacts is not None
             and contacts.soft_contact_max > 0