[geometry] Add optional MuJoCo multi-point contact for convex meshes

MODULE.bazel

@@ -328,6 +328,7 @@ use_repo(
     "libzip_internal",
     "metis_internal",
     "msgpack_internal",
+    "mujoco_internal",
     "nanoflann_internal",
     "nlopt_internal",
     "onetbb_internal",

bindings/generated_docstrings/geometry.h

@@ -990,6 +990,31 @@ Collision Detection.
     Currently margin only applies to *compliant* hydroelastic contact
     and it does not affect point contact.)""";
         } margin;
+        // Symbol: drake::geometry::DefaultProximityProperties::point_contact_algorithm
+        struct /* point_contact_algorithm */ {
+          // Source: drake/geometry/scene_graph_config.h
+          const char* doc =
+R"""((Experimental) Selects the narrowphase algorithm used to compute point
+contact for this geometry. There are two valid options: -
+"single_point": the existing behavior; every penetrating geometry pair
+reports exactly one contact point (the deepest one). -
+"mujoco_multipoint": pairs where *both* geometries are Convex or Mesh
+shapes (and both request this algorithm) report a contact manifold of
+up to four contact points, computed with MuJoCo's native convex
+collision detection (GJK/EPA plus contact-polygon clipping over the
+meshes' convex hulls); each point carries its own depth. Face-face and
+edge-face contacts gain the extra points; vertex-face contacts still
+report one point. All other geometry pairs keep the "single_point"
+behavior. A mesh whose convex hull has no volume (for example, a
+planar mesh) reports a single point even when it selects this
+algorithm.
+
+This affects QueryObject∷ComputePointPairPenetration() and the point
+contact fallback of QueryObject∷ComputeContactSurfacesWithFallback()
+for double-valued scene graphs. Scene graphs of other scalar types
+ignore the property; those scalar types do not support penetration
+queries between Convex or Mesh shapes in the first place.)""";
+        } point_contact_algorithm;
         // Symbol: drake::geometry::DefaultProximityProperties::point_stiffness
         struct /* point_stiffness */ {
           // Source: drake/geometry/scene_graph_config.h
@@ -1049,6 +1074,7 @@ R"""(See also:
             std::make_pair("hunt_crossley_dissipation", hunt_crossley_dissipation.doc),
             std::make_pair("hydroelastic_modulus", hydroelastic_modulus.doc),
             std::make_pair("margin", margin.doc),
+            std::make_pair("point_contact_algorithm", point_contact_algorithm.doc),
             std::make_pair("point_stiffness", point_stiffness.doc),
             std::make_pair("relaxation_time", relaxation_time.doc),
             std::make_pair("resolution_hint", resolution_hint.doc),
@@ -6129,6 +6155,16 @@ values reported here are confirmed, observed worst case answers.
     For Mesh shapes, their convex hulls are used in this query. It is
     not* computationally efficient or particularly accurate.
 
+Note:
+    By default every penetrating geometry pair reports exactly one
+    point pair (its deepest point). A Convex-Convex, Convex-%Mesh, or
+    Mesh-%Mesh pair where *both* geometries select the (experimental)
+    "mujoco_multipoint" point contact algorithm (see
+    DefaultProximityProperties∷point_contact_algorithm) instead
+    reports a contact manifold of *up to four* point pairs for ``T`` =
+    ``double``, each with its own depth, ordered deterministically for
+    fixed poses.
+
 Raises:
     RuntimeError if a Shape-Shape pair is in collision and indicated
     as ``throws`` in the support table above.)""";

bindings/generated_docstrings/geometry_proximity.h

@@ -49,6 +49,8 @@
 // #include "drake/geometry/proximity/mesh_to_vtk.h"
 // #include "drake/geometry/proximity/mesh_traits.h"
 // #include "drake/geometry/proximity/meshing_utilities.h"
+// #include "drake/geometry/proximity/mujoco_ccd_mesh_data.h"
+// #include "drake/geometry/proximity/mujoco_ccd_penetration.h"
 // #include "drake/geometry/proximity/obb.h"
 // #include "drake/geometry/proximity/obj_to_surface_mesh.h"
 // #include "drake/geometry/proximity/plane.h"

bindings/pydrake/geometry/test/scene_graph_test.py

@@ -485,6 +485,7 @@ def test_scene_graph_config(self):
             hunt_crossley_dissipation=7,
             relaxation_time=None,  # Test optionality.
             point_stiffness=9,
+            point_contact_algorithm="mujoco_multipoint",
         )
         param_init_scene_graph = mut.SceneGraphConfig(
             default_proximity_properties=param_init_props
@@ -493,6 +494,7 @@ def test_scene_graph_config(self):
         got_props = param_init_scene_graph.default_proximity_properties
         self.assertEqual(got_props.relaxation_time, None)
         self.assertEqual(got_props.point_stiffness, 9)
+        self.assertEqual(got_props.point_contact_algorithm, "mujoco_multipoint")
 
     @numpy_compare.check_all_types
     def test_scene_graph_renderer_with_context(self, T):

examples/multibody/multipoint_contact/BUILD.bazel

@@ -0,0 +1,15 @@
+load("//tools/lint:lint.bzl", "add_lint_tests")
+load("//tools/skylark:drake_py.bzl", "drake_py_binary")
+
+drake_py_binary(
+    name = "falling_box",
+    srcs = ["falling_box.py"],
+    deps = ["//bindings/pydrake"],
+    add_test_rule = True,
+    test_rule_args = [
+        "--target_realtime_rate=0",
+        "--simulation_time=0.5",
+    ],
+)
+
+add_lint_tests()

examples/multibody/multipoint_contact/README.md

@@ -0,0 +1,107 @@
+# Multi-point contact between convex shapes
+
+This example shows what the `"mujoco_multipoint"` point contact algorithm
+changes, by dropping two boxes onto a slab side by side: the box on the left
+touches the slab through Drake's default single contact point, the box on the
+right through a four-point contact manifold.
+
+## Background
+
+Drake's point contact model resolves each colliding geometry pair to a single
+contact point, the point of deepest penetration, even when two flat faces
+overlap. A box resting on a table is then supported at one point and cannot
+be in static equilibrium there: the support point wanders across the bottom
+face and the box rocks. The usual workaround is to decorate a body with
+several small spheres so that it touches through several points (see
+`examples/multibody/cylinder_with_multicontact`).
+
+The proximity property `("material", "point_contact_algorithm")` selects an
+alternative narrowphase for point contact. Its value `"mujoco_multipoint"`
+uses MuJoCo's native convex collision detection (GJK/EPA followed by clipping
+the two touching faces against each other) to report a contact manifold of up
+to four points, each with its own penetration depth, for pairs of `Convex` or
+`Mesh` shapes. Face-face and edge-face contacts gain the extra points;
+vertex-face contacts still report a single point. Two rules matter for setting
+up a scene:
+
+- *Both* geometries of a pair must select `"mujoco_multipoint"`; a pair with
+  one geometry left at the default `"single_point"` reports one point.
+- Only `Convex` and `Mesh` shapes are eligible. A `Box` shape is not, so this
+  example models its boxes as `Convex` shapes built from an in-memory OBJ file
+  listing the eight corners (`make_convex_box()` in `falling_box.py`).
+
+The property can be set per geometry, as this example does, or for a whole
+scene at once through
+`SceneGraphConfig.default_proximity_properties.point_contact_algorithm`.
+
+## The scene
+
+- A 1 m x 0.5 m x 0.05 m slab, a box modeled as a `Convex` shape, welded to
+  the world with its top face at z = 0. It selects `"mujoco_multipoint"`.
+- Two identical 20 cm x 20 cm x 4 cm boxes of 0.5 kg, also `Convex` shapes.
+  The orange box on the left selects `"single_point"`; the blue box on the
+  right selects `"mujoco_multipoint"`.
+- Both boxes start level, 5 cm above the slab, and fall onto it.
+
+The plant uses `contact_model = "point"` with a 1 ms discrete time step.
+
+## Running
+
+```
+bazel run //examples/multibody/multipoint_contact:falling_box
+```
+
+Open the Meshcat URL printed at startup. The contact forces are drawn as
+arrows, one per contact point, so the left box shows a single arrow and the
+right box shows four once it has landed. The console prints, every half
+second, the number of contact points each box has with the slab, the height
+of the box center above the slab in millimeters (20.00 for a box at rest on
+the slab), and the box's angular speed in rad/s:
+
+```
+               single_point box               multipoint box
+   time   contacts   height    |w|     contacts   height    |w|
+   0.50          1    17.46   0.242          4    20.00   0.000
+   1.00          1    16.31   0.237          4    20.00   0.000
+   ...
+   5.00          1    14.50   0.228          4    20.00   0.000
+```
+
+Use `--simulation_time`, `--target_realtime_rate`, `--time_step`, and
+`--drop_height` to change the run.
+
+## What to look for
+
+Both boxes reach the slab about a tenth of a second after the start. From
+then on the two halves of the table diverge.
+
+The single-point box reports one contact, and that contact hops between the
+corners of its bottom face from one time step to the next. The box never comes
+to rest: it rocks at about 0.23 rad/s for the whole run, creeps sideways and
+in yaw by a few degrees over several seconds, and sinks into the slab. Its
+center, which would sit 20.00 mm above the slab at rest, reads 16.3 mm after
+one second and 14.5 mm after five, so the box has sunk 5.5 mm into a slab it
+should be resting on. (The contact stiffness alone would account for a few
+micrometers of penetration under the box's 5 N weight.) In Meshcat the single
+force arrow jumps around the bottom face.
+
+The multipoint box reports four contacts, one at each corner of its bottom
+face, and rests at exactly 20.00 mm with zero angular velocity from the moment
+it lands. Meshcat shows four steady arrows.
+
+The size of the single-point artifact depends on the time step; the manifold
+result does not. With `--time_step=0.01` the single-point box falls straight
+through the slab within half a second while the multipoint box still rests on
+it. With `--time_step=0.0001` the single-point box settles to within 0.01 mm of
+its rest height but still wobbles at 0.02 to 0.04 rad/s.
+
+Why one point is not enough: a level box supported at a single point is in
+equilibrium only if that point lies directly under its center of mass. The
+deepest point of a face-face overlap is a corner (whichever one a tiny
+numerical tilt favors), so the normal force there tips the box toward the
+opposite corner, which then becomes the deepest point, and the box rocks from
+corner to corner. Each exchange of support corners lets the center of mass
+fall a little before the contact catches it, so the box ratchets downward,
+and a larger time step makes each exchange coarser and the sinking faster.
+Four contact points spanning the overlap region support the box the way a
+table does.