Feature/surface velocity

[juanmed]

Background

This is the first PR to address #19599. This feature consists on modelling a conveyor belt by defining a velocity at the surface of an each we want to show such conveyor belt effect. Most of the design considerations where discussed in the issue mentioned above, here is a summary of those considerations:

• The surface velocity through three quantities:

  • Velocity normal n_v : a 3-vector expressed in the frame of the geometry we want to have a surface velocity
  • Surface normal n: a 3-vector representing the normal to the surface at some point on the geometry
  • Surface speed s : the magnitude of the surface velocity
    The surface velocity v_s vector is orthogonal to both n_v and n and computed with: v_s = s * n_v x n. (x is the cross product)

• There are 3 ways to supply the quantities required to define a surface velocity for a geometry. These are the velocity normal and surface speed (surface normal depends on the Shape of the geometry):

  • Input ports: created on demand when calling a member function like MultibodyPlant::DeclareSurfaceVelocityInput(const GeometryId id, const Vector3<T>& default_velocity_normal, const T& default_belt_speed);. The value of the velocity normal and surface speed can change dynamically with the simulation.
  • Default values: passed in through the corresponding parameters in the function above
  • SDF/URDF: define both as ProximityParameters which are read from drake:proximity_parameters tag in sdf/urdf files.

• Surface velocity is defined per-geometry and not per-body. This is different from the original discussion. The reason behind this decision is that the definition of the surface velocity through SDF/URDF ProximityParameters is done per-geometry. So it makes sense to keep it consistent across APIs. Otherwise, when using SDF/URDF the user would work on a per-geometry basis, and have to think on a per-body basis when using other methods.

Description

This PR implements the basic functionality to define a surface velocity. Namely:

• Define and load drake:surface_speed and drake:surface_velocity_normal parameters as ProximityParameters. See proximity_properties.h/cc,
• Load those parameters from urdf or sdf files. Refer to detail_common.c, detail_sdf_parser.cc and detail_urd_parser.cc
• Functions to compute the surface normal at some point over the surface defined by the Shape of some geometry defined by GeometryId. shape_specification.cc/h
• Compute surface velocity for point-contact and the continuous (multibody_plant.cc/h) and discrete (discrete_update_manager.h) plants

Other changes:

• Breakup some dependencies on shape_specification.h to avoid circular dependencies. In particular for bvh, bv targets in geometry folder. This allows using distance computation functions within Convex and Mesh shapes.

The math behind modelling the surface velocity is mostly based on [1], [2], and the specific integration of surface velocity is described in the following slides:

[1] A. Castro, F. Permenter, and X. Han, “An Unconstrained Convex Formulation of Compliant Contact,” June 23, 2022, arXiv: arXiv:2110.10107. doi: 10.48550/arXiv.2110.10107.
[2] A. Castro, X. Han, and J. Masterjohn, “Irrotational Contact Fields,” July 15, 2025, arXiv: arXiv:2312.03908. doi: 10.48550/arXiv.2312.03908.

Tests

Example using the features introduced in this PR bazel run //examples/conveyor_belt:conveyor_belt

Unit tests:
bazel run //geometry:shape_specification_test : verify computation of normal at a point on the surface of a geometry::Shape
bazel run //multibody/parsing:detail_sdf_geometry_test and bazel run //multibody/parsing:detail_sdf_geometry_test : test reading proximity parameters from urdf and sdf files.
bazel run //multibody/plant:surface_velocity_test: check basic surface velocity features

Still pending
Next PR covers the following:

• Compute surface velocity for hydroelastic-contact
• Input ports to change surface speed and velocity normal dynamically

I can merge that work here if reviewers would like to see the whole work in one go. However the size of the PR will increase significantly.

---

This change is 

[SeanCurtis-TRI]

multibody/plant/multibody_plant.cc line 2362 at r1 (raw file):

    const Vector3<T> v_WBc =
        vc.get_V_WB(bodyB_mobod_index).Shift(-p_CoBo_W).translational();
    const Vector3<T> v_AcBc_W = v_WBc - v_WAc;

This is the relative velocity of the two bodies at the point of contact. Your surface speeds should contribute to this vector. Currently, this only includes the velocities of the bodies.

Remember that the surface speed (and direction of travel) is defined in the body's frame. So, when combining it here you'll need to map it to the body frame so that all of the vectors meaningfully combine.

[jwnimmer-tri]

I'll pick +@amcastro-tri to help shepherd this along, but of course feel free to delegate it, too.

[juanmed]

@jwnimmer-tri Thanks for taking a look here, any advice is appreciated.

At this point the next thing for me is looking at how to integrate this with the SAP solver as suggested by @amcastro-tri to verify simulation is sensible. I am not familiar with SAP so any indications would help me get started faster. Aftewards I would look into creating the input port to pass in the surface velocity and finally integrating the surface velocity in the hydroelastic model calculations. Implementation-wise any feedback would be great as I definitely do not have the bigger picture of where to implement what... so far I have gone with my intuition after reading relevant sections.

What I have so far is:

• Pass in surface speed and velocity plane normal as proximity parameters in sdf file (still to look at urdf)
• Read those proximity parameters in MultibodyPlant and compute the surface velocity vector
• Correctly combine the surface velocity with separation velocity as suggested above by @SeanCurtis-TRI
• WIP unit testing these features.

I am also leaving concrete comments on different sections.

[jwnimmer-tri]

FYI for posting comments in the code, please open the "Reviewable" link that got automatically added to the PR overview and use that tool for discussion. In Drake, we don't use GitHub's (impoverished) code review tool.

[juanmed]

Ok. Will use that, thanks for the note.

[juanmed]

multibody/plant/multibody_plant.cc line 2362 at r1 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

This is the relative velocity of the two bodies at the point of contact. Your surface speeds should contribute to this vector. Currently, this only includes the velocities of the bodies.

Remember that the surface speed (and direction of travel) is defined in the body's frame. So, when combining it here you'll need to map it to the body frame so that all of the vectors meaningfully combine.

Done.

[juanmed]

Seems to be working now with Box and continuous plant.

conveyor_belt-2025-07-23_13.54.26.mp4

But there some weird artifacts still, like this box that lands at an angle on the conveyor belt: at first moves in the correct direction but later moves backward and falls through the conveyor (?). Not sure why this happens yet.

conveyor_belt-2025-07-23_14.09.49.mp4

[amcastro-tri]

Nice!
I caught the glimpse of a green force arrow for a split second. That means you have "point contact". Point contact for boxes is awful and can lead to this kind of artifacts. I recommend you switch to using the hyudroelastic contact model.

Reviewable status: 12 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @juanmed)

[juanmed]

A bit more progress: for a discrete plant and using SAP solver.

conveyor_belt-2025-07-25_00.13.16.mp4

The behavior when the box is initially at an angle is visually better when compared to the continuous plant:

conveyor_belt-2025-07-25_00.21.37.mp4

[amcastro-tri]

Excellent progress! I love this! I placed a few comments below. Aftear reading your apporach (very clever!) I realize that probably there's a yet even simpler solution (modify v_AcBc_C to incoroporate surface velocities). PTAL

Reviewed 2 of 16 files at r6.
Reviewable status: 16 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @juanmed)

---

multibody/contact_solvers/contact_configuration.h line 87 at r6 (raw file):

  // it should be already added to vn on construction of this
  // ContactConfiguration.
  Vector3<T> vt_b{0., 0., 0.};

Correct me if I am wrong, we could probably describe this in physical terms more easily.
I believe this is the relative velocity of the belt with respet to the body they belt is attached to, expressed in the contact frame.
Is that correct?

---

multibody/plant/multibody_plant.cc line 1713 at r6 (raw file):

  } else {
    return std::nullopt;

question in #19599 it was suggested to use input ports instead. That'd allow chaning the belt's velocity (say on/off) from an external system. If placed in a proximity property, they you'd need to go outside the main simulation loop, modify the proximity properties, and then start the simulation again.

Also, it turns out that retrieving properties from ProximityProperties is extremely inneficient (historic reasons). Quering a port is a lot faster.

---

multibody/contact_solvers/sap/sap_hunt_crossley_constraint.cc line 170 at r6 (raw file):

  // Computations dependent on vc.
  data.vc = vc + configuration_.vt_b;
  data.vn = data.vc.z();

nice! do you foressee any case in which vt_b.z() is different form zero? adding a normal velocity could cause very strange behavior.

---

multibody/plant/discrete_update_manager.cc line 723 at r6 (raw file):

    // Contact velocity stored in the current context (previous time step).
    const Vector3<T> v_AcBc_W = Jv_AcBc_W * v;
    const Vector3<T> v_AcBc_C = (R_WC.transpose() * v_AcBc_W);

your implementation made me think. I believe we don't need to modify the SAP constraints at all, all we need to do is to modify the defintion of the slip velocity vt to include the surface velocities. That is, what we want is to impose the slip velocity to be zero, subject to Coulomb's law. When the object has a surfacce velocity (say a belt wrapping it), we no longer care about the velocity of the underlying object, only the velocity of the belt. That is, I belive all we need to add is to incorporate the surfaces velocities (as you do here) "directly" into the definition of v_AcBc_W. Then the constraint do not even need to be modified.

We would also need to do this for hydroelastic contact in "AppendDiscreteContactPairsForHydroelasticContact()".

Even more, we could have a single entry point, say DiscreteUpdateManager::AddSurfaceVelocities(const Context& context, std::vector<DiscreteContactPair>* pairs) that takes processes all discrete pairs in a single place to add the contact velocities.

[juanmed]

multibody/contact_solvers/contact_configuration.h line 87 at r6 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

Correct me if I am wrong, we could probably describe this in physical terms more easily.
I believe this is the relative velocity of the belt with respet to the body they belt is attached to, expressed in the contact frame.
Is that correct?

Not really, this is relative surface between the two bodies at the point of contact: v_Bc_surf - v_Ac_surf, expressed in the contact frame. In that sense it is not necessarily any of the surface velocities of the bodies but the vector sum of both (recall both bodies could have each a conveyor belt).

[juanmed]

multibody/contact_solvers/contact_configuration.h line 87 at r6 (raw file):

Previously, juanmed (Juan) wrote…

Not really, this is relative surface between the two bodies at the point of contact: v_Bc_surf - v_Ac_surf, expressed in the contact frame. In that sense it is not necessarily any of the surface velocities of the bodies but the vector sum of both (recall both bodies could have each a conveyor belt).

I tried to improve the wording here following your description a bit.

[juanmed]

multibody/contact_solvers/sap/sap_hunt_crossley_constraint.cc line 170 at r6 (raw file):
In general, vt_b.z() is non-zero. The only case where we could assume it is zero is in the case where the surface normal and the contact normal are the same: the surface velocity does not have a but in general this is not the case. Recall that in general the surface normal and contact normal need to be the same.

adding a normal velocity could cause very strange behavior.

In the experiments I have made so far, I did not see any strange behavior. If you have any tests in mind I am happy to try them.

[juanmed]

multibody/plant/discrete_update_manager.cc line 723 at r6 (raw file):
I agree with your ideas and initially I was looking for a way to not modify the constraints. But realized v_AcBc_W is never passed downstream (only its normal component vn0), and figured there is no other choice but passing v_AcBc_C_ss as a new member of ContactConfiguration. Correct me if I am wrong, maybe I missed the mechanism by which v_AcBc_W is passed down.

When the object has a surfacce velocity (say a belt wrapping it), we no longer care about the velocity of the underlying object, only the velocity of the belt.

Maybe I misunderstood but I think this is incorrect: the velocity of the object continues to be relevant even when there is a surface velocity defined. Imagine an articulated robot arm with a small conveyor belt end-effector attached to it. When the conveyor belt interacts with any object, it does so at a velocity that is the vector sum of the velocity of the end-effector link and the surface velocity of the conveyor belt at the point of contact.

[juanmed]

multibody/plant/discrete_update_manager.cc line 723 at r6 (raw file):

Even more, we could have a single entry point, say DiscreteUpdateManager::AddSurfaceVelocities(const Context& context, std::vector<DiscreteContactPair>* pairs) that takes processes all discrete pairs in a single place to add the contact velocities.

This sounds good. The current PR will limit to point-contact so will look at this idea in the next PR.

[juanmed]

multibody/plant/multibody_plant.cc line 1713 at r6 (raw file):

Previously, amcastro-tri (Alejandro Castro) wrote…

question in #19599 it was suggested to use input ports instead. That'd allow chaning the belt's velocity (say on/off) from an external system. If placed in a proximity property, they you'd need to go outside the main simulation loop, modify the proximity properties, and then start the simulation again.

Also, it turns out that retrieving properties from ProximityProperties is extremely inneficient (historic reasons). Quering a port is a lot faster.

Yes, an input port is definitely to be implemented in a follow up PR (I already have this working). For this PR I will limit to using Proximity Properties because:

• Using drake proximity properties allows to define a surface velocity in the .sdf / .urdf , making it easier to use this feature for a simple conveyor without caring about defining and connecting ports.
• The input port, if exists, will always override any values passed by ProximityProperties

This is what I understood from the whole conversation regarding interfacing to this feature: #19599 (comment)

[juanmed]

Thanks for taking a look. I replied in each of your comments.

[juanmed]

@amcastro-tri This is ready for review. \

Some linting tests are failing in CI, possibly because my tools are not configured as drake expects. Is there any developer setup file I can follow to format/lint this as expected?

[rpoyner-tri]

I pushed some format/build fixups. It appears there are still other problems.

Reviewable status: 6 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @juanmed)

[juanmed]

I pushed some format/build fixups. It appears there are still other problems.

Reviewable status: 6 unresolved discussions, LGTM missing from assignee amcastro-tri, needs platform reviewer assigned, needs at least two assigned reviewers, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @juanmed)

Thank you! Will check the errors.

[juanmed]

Bump

[juanmed]

Bump

[juanmed]

Bump

[jwnimmer-tri]

One or both of @SeanCurtis-TRI and @joemasterjohn will help out here.

(Also please merge up to latest master, resolving conflicts, and push an update.)

[SeanCurtis-TRI]

@juanmed Apologies, I just realized you're waiting on us. I'll go over this and get back to you in the next day or so. Thanks for your patience.

[SeanCurtis-TRI]

Sorry for how longs this has taken. I've had to digest everything here anew and there's a fair amount going on. I'm not finished all the way through but I feel I've hit something big enough that it would be helpful to get that back to you now. See the next discussion about what normal to use in the contact calculations.

I'll continue through and publish the remaining comments tomorrow.

@SeanCurtis-TRI reviewed 15 files and made 23 comments.
Reviewable status: 25 unresolved discussions, LGTM missing from assignees SeanCurtis-TRI(platform),joemasterjohn, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on juanmed).

---

multibody/contact_solvers/contact_configuration.h line 87 at r6 (raw file):

Previously, juanmed (Juan) wrote…

I tried to improve the wording here following your description a bit.

I certainly by the argument that this bias term is not necessarily the relative velocity w.r.t. the body.

That said, I'd recommend the following:

1. Call it Vector3<T> vt_b as a bias for the tangent component of the contact velocity, expressed in the contact frame. So, vt_b.z() is always zero.
2. In the documentation, give the surface_velocity example as to why this bias term might not be zero (maybe in the future, we'll have other reasons for non-zero bias, but, for now, it's worth giving the reader a hint as to why this exists.

---

multibody/contact_solvers/sap/sap_hunt_crossley_constraint.cc line 170 at r6 (raw file):

Previously, juanmed (Juan) wrote…

In general, vt_b.z() is non-zero. The only case where we could assume it is zero is in the case where the surface normal and the contact normal are the same: the surface velocity does not have a but in general this is not the case. Recall that in general the surface normal and contact normal need to be the same.

adding a normal velocity could cause very strange behavior.

In the experiments I have made so far, I did not see any strange behavior. If you have any tests in mind I am happy to try them.

Based on all previous recommendations, this new code is fine, because the newly named vt_b would have a zero in the third element. But, the key is, it would always be zero.

---

multibody/plant/discrete_update_manager.cc line 723 at r6 (raw file):

When the object has a surfacce velocity (say a belt wrapping it), we no longer care about the velocity of the underlying object, only the velocity of the belt.

If I can give my own interpretation to this. I don't believe he means that the velocity of the body doesn't matter. But it only matters in so far as it changes the velocity of the point on the belt w.r.t. an inertial frame. Once we have that quantity, the body velocities hold no further significance.

At the end of the day, what we care about is the velocity at the surface which is a combination of surface velocity and body velocity.

---

a discussion (no related file):
I have a fundamental question.

You've characterized the "belt velocity" (if I can use that term) as a function of point on surface, normal at point, and velocity normal.

What I'm confused about is the selection of the point and surface normal. I'll claim that (for contact), the point and normal should always simply be the contact point and the contact normal, respectively.

Arguments

Physically, objects don't intersect. When real objects make contact, the contact point(s) lie(s) on the surface (even if the object deforms to achieve this). The contact normal is, by definition, the surface normal at that contact. So, in the real world, there is no distinction.

We've got compliant models of contact. Generally, we know the contact point will not lie on the surface of either of the intersecting bodies. But, what our models do, in principle, is approximate the minute deformations of the objects (increasing the elastic potential energy in the materials) and impart a force. In that approximation, the contact point truly lies on that approximate deformed object's surface and the contact normal is simultaneously perpendicular to both object's deformed surfaces (such as they are).

Furthermore, projecting a contact point that lies inside a shape onto its surface can produce horribly discontinuous results. It's most obvious in the case of boxes. If a point lies near the boundary of the voronoi regions on the interior of a box, a small perturbation can lead the "nearest" point to jump from one face to another. That is definitely not desirable. We want toe contact data to be smooth w.r.t. configuration.

The proposed conveyor belt model is not truly a surface velocity, but a material velocity field. Sampled on the surface of a geometry, it is properly the surface velocity. However, we still define velocities away from the surface that converge to the surface velocity as the contact gets shallower and shallower.

Implications

1. First, there will never be any component of the conveyor belt velocity that is in the contact normal direction. That is good and proper. The conveyor belt should never attract or repulse and object.
2. It reduces our bookkeeping and implementation. We have few point and normals to account for (or generate).

---

geometry/proximity_properties.cc line 16 at r8 (raw file):

const char* const kPointStiffness = "point_contact_stiffness";
const char* const kSurfaceVelocityGroup = "surface_velocity";
const char* const kSurfaceSpeed = "surface_speed";

nit GIven you already have a group named "surface_velocity", adding the name "surface" to the property name is redundant. Think about the property as the tuple (group, name) and contrast:

(surface_velocity, speed)
(surface_velocity, normal)

with

(surface_velocity, surface_speed)
(surface_velocity, surface_normal)

---

geometry/shape_specification.h line 0 at r8 (raw file):
Generally, shape_specificaiton is just that specification. It's generally tried to avoid calculations.

That said, based on my assertion that the belt velocity should be based strictly on contact point and normal, all of this goes away.

---

geometry/shape_specification.h line 826 at r8 (raw file):

// returns no value.
template <typename T>
std::optional<Eigen::Vector3<T>> GetNormalAtPoint(const Shape& shape,

nit: This API is especially suspect. As near as I can tell, it's only called in MultibodyPlant. At that point, what is being passed is the contact point in the geometry frame. The contact point will generally not be on the surface which means this function would be obliged to return nullopt, generllly (which in tern means returning a zero surface velocity).

But, again, this is just another point in support of, "use the contact normal to determine belt velocity".

---

geometry/proximity/aabb.h line 0 at r8 (raw file):
This feels like a merging problem. This should be part of master and not part of this PR. Try merging from master again and we'll see what disappears.

---

geometry/proximity/test/hydroelastic_internal_test.cc line 18 at r8 (raw file):

#include "drake/geometry/proximity/make_sphere_field.h"
#include "drake/geometry/proximity/make_sphere_mesh.h"
#include "drake/geometry/proximity/proximity_shape_utilities.h"

nit: Unused by this PR.

---

geometry/proximity/test/make_capsule_mesh_test.cc line 6 at r8 (raw file):

#include "drake/common/eigen_types.h"
#include "drake/geometry/proximity/proximity_shape_utilities.h"

nit: Unused by this PR.

---

geometry/proximity/test/mesh_half_space_intersection_test.cc line 17 at r8 (raw file):

#include "drake/common/test_utilities/expect_throws_message.h"
#include "drake/geometry/proximity/contact_surface_utility.h"
#include "drake/geometry/shape_specification.h"

nit: Unused by this PR.

---

geometry/proximity/test/mesh_intersection_test.cc line 16 at r8 (raw file):

#include "drake/common/test_utilities/expect_throws_message.h"
#include "drake/geometry/geometry_ids.h"
#include "drake/geometry/shape_specification.h"

nit: Unused by this PR.

---

geometry/proximity/test/mesh_plane_intersection_test.cc line 18 at r8 (raw file):

#include "drake/common/unused.h"
#include "drake/geometry/proximity/contact_surface_utility.h"
#include "drake/geometry/shape_specification.h"

nit: Unused by this PR.

---

geometry/proximity/test/proximity_utilities_test.cc line 14 at r8 (raw file):

namespace {

using Eigen::Vector3d;

nit: Unused by this PR.

---

geometry/test/shape_to_point_cloud.h line 0 at r8 (raw file):
Does this belong in this PR?

---

multibody/contact_solvers/contact_configuration.h line 86 at r8 (raw file):

  // where vc is the contact velocity, J the contact jacobian and v the vector
  // of generalized velocities. In practice, is used to model a velocity at the
  // contact point, such as when an imaginery conveyor belt is wrapped around

nit typo

Suggestion:

imaginary

---

multibody/contact_solvers/sap/sap_friction_cone_constraint.cc line 101 at r8 (raw file):

  data.mutable_vc() = vc + configuration_.v_b;
  data.mutable_y() =
      data.R_inv().asDiagonal() * (data.v_hat() - vc - configuration_.v_b);

nit: this would be more robust and readable with a slight change.

Suggestion:

  data.mutable_y() =
      data.R_inv().asDiagonal() * (data.v_hat() - data.vc());

---

multibody/contact_solvers/sap/sap_hunt_crossley_constraint.cc line 57 at r8 (raw file):

  p.dt = time_step;
  const T& fe0 = configuration_.fe;
  const T& vn0 = configuration_.vn + configuration_.v_b.z();

Based on all previous recommendations, there would be no bias term in the normal direction.

---

multibody/parsing/detail_common.cc line 0 at r8 (raw file):
BTW I'm skpping parsing for now. I'll revisit it after we've cleaned up the math implementation.

BTW Generally, it's easier to review these kinds of PRs if you decompose it a bit more. One PR that gets in the functionality, one that handles the parsing. It keeps the PR easier to consume.

---

multibody/plant/deformable_driver.cc line 588 at r8 (raw file):

          .phi0 = phi0,
          .vn0 = v_AcBc_Cz,
          .v_b = {},

nit: This should be Vector3<T>::Zero(). Default initialized would leave garbage in this field.

---

multibody/plant/discrete_contact_pair.h line 96 at r8 (raw file):

  T vn0{0.0};
  /* Tangential velocity bias term, used to model a velocity at the surface*/
  Vector3<T> v_b;

nit: Same not here as for the sap-based contact data.

---

multibody/plant/discrete_update_manager.cc line 709 at r8 (raw file):

    // Get surface velocity at Ca relative to A in coordinates of A and
    // transform to world frame W.

BTW There's a canonical Drake way to say this.

(Applies, to the next as well.)

Suggestion:

    // Get surface velocity at Ca relative to Ao measured in A and expressed in
    // World.

[SeanCurtis-TRI]

This badly needs a remerge with master:

1. The build infrastructure has moved on (new version of bazel).
   • Without doing this, it's incredibly cumbersome to build this locally if I have otherwise moved my local drake branch to master and updated my build infrastructure.
2. There are now some heavy conflicts between some of this PR and master: the changes to shape_specification (such that it depends on downstream elements like bvh/obb). Add Aabb-{plane|half space} intersection and missing bindings #24001 introduced conflicting implementations with what has been added in this PR. That implementation explicitly depends on shape specification (by naming HalfSpace).

I have a couple of ideas for testing my assertions about the potential defects in implementation but, in order to verify them, I need to be able to build. So, I'm going to merge master and reconcile conflicts without changing the spirit of the implementation despite the fact that I believe we'll be removing most (if not all) of the geometry changes subsequently.

[SeanCurtis-TRI]

@juanmed I've more or less finished the full review of this (well, as full as I'm going to do without the discussion below about where your surface normal comes from getting resolved). There is still a significant model problem (see discrete_update_manager.cc)

Here's the big question: you were industriously working on this last fall before we disappeared. I don't know what you're currently up to. If you're still interested in landing this in Drake (with the attendant follow up of landing including hydroelastic contact) then I'm happy to help you push this over the finish line.

If, on the other hand, you've moved on, then I'll consider what we'd like to do with this functionality. Either way, if you could let us know what you're up for, we'd appreciate. Again, apologies for leaving you hanging.

@SeanCurtis-TRI reviewed 4 files and made 10 comments.
Reviewable status: 35 unresolved discussions, LGTM missing from assignees SeanCurtis-TRI(platform),joemasterjohn, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on juanmed).

---

multibody/plant/discrete_update_manager.cc line 703 at r8 (raw file):

        math::RotationMatrix<T>::MakeFromOneVector(nhat_AB_W, 2);

    const RigidTransform<T>& X_WA =

nit: In the PR description, you've claimed that the surface velocity normal is expressed in the geometry frame. So, for body frame A, the frame of the geometry G affixed to it would have the pose X_AG. Here, we're computing X_WA. X_WA is only equal to X_WG if X_AG = I. That is not generally the case.

You can see this easily if you move the pose in conveyor_belt.sdf from the <link> level to the <collision> and <visual> level. The scene will look the same and the box and conveyor belt will collide in the same configuration, but there will be no conveyor belt action.

---

multibody/plant/discrete_update_manager.cc line 719 at r8 (raw file):

        plant().GetSurfaceVelocity(pair.id_B, inspector, X_WB, pair.p_WCb);
    // Relative separation velocity due to surface velocity in contact frame C.
    const Vector3<T> v_AcBc_C_ss = R_WC.transpose() * (v_BCb_W_ss - v_ACa_W_ss);

This expression is hugely suspicious.

As reasoned about here, you have two different velocities measured in different frames (although they are expressed in a common frame). Velocities measured in different frames simply can't be meaningfully subtracted.

This is like subtracting the velocity of a person walking down a train aisle from the velocity of a person walking down a plane aisle. If you want to talk about the relative velocities of the two people, relative to the ground, you must account for the velocities of the plane and train w.r.t. the ground. This shortcuts straight across that (and I feel contributes to issues pointed out elsewhere).

I've got a geometric way to reason about this that makes things much simpler. I've tested it against the various problem cases I've called out above and it seems to do what is expected.

I'll wait to see what we decide to do about this PR.

---

examples/conveyor_belt/conveyor_belt.sdf line 99 at r8 (raw file):

        <drake:proximity_properties>
          <drake:mu_dynamic>0.7</drake:mu_dynamic>
          <drake:mu_static>0.7</drake:mu_static>

BTW You've applied moving surfaces to all the shapes except this box. The others behave as expected, they slide across the ground.

I should be able to apply velocity to this box as well. The interaction between this box and the conveyor belt it falls on should combine. When I copy the conveyor belt's surface velocity configuration onto this box, I don't get additive sliding effects -- the box gets sucked through the table.

Again, this points to a couple of things:

1. The surface velocity should always be tangent to the contact normal. It should only contribute to friction calculation and never to normal forces.
2. The calculation of relative velocity is probably adding to the problems here.

---

examples/conveyor_belt/BUILD.bazel line 0 at r8 (raw file):
If we keep the example, we should have a Readme.md file in here explaining what is being exampled.

---

examples/conveyor_belt/BUILD.bazel line 21 at r8 (raw file):

        "//systems/analysis:simulator",
        "//systems/framework:diagram",
        "//systems/framework:leaf_system",

To support alternate visualization

Suggestion:

        "//geometry:scene_graph",
        "//multibody/parsing",
        "//multibody/plant",
        "//systems/analysis:simulator",
        "//systems/framework:diagram",
        "//systems/framework:leaf_system",
        "//visualization",

---

examples/conveyor_belt/conveyor_belt.cc line 10 at r8 (raw file):

#include "drake/multibody/plant/multibody_plant.h"
#include "drake/systems/analysis/simulator.h"
#include "drake/systems/framework/diagram_builder.h"

To change visualization

Suggestion:

#include "drake/geometry/meshcat.h"
#include "drake/geometry/scene_graph.h"
#include "drake/multibody/parsing/parser.h"
#include "drake/multibody/plant/multibody_plant.h"
#include "drake/systems/analysis/simulator.h"
#include "drake/systems/framework/diagram_builder.h"
#include "drake/visualization/visualization_config.h"
#include "drake/visualization/visualization_config_functions.h"

---

examples/conveyor_belt/conveyor_belt.cc line 39 at r8 (raw file):

  cparams.newtons_per_meter = 60.0;
  multibody::meshcat::ContactVisualizerd::AddToBuilder(&builder, plant, meshcat,
                                                       std::move(cparams));

nit: The superior way to add visualization is simpler:

See the other comments for the requisite changes.

Suggestion:

  visualization::ApplyVisualizationConfig(
      visualization::VisualizationConfig{
          .default_proximity_color = geometry::Rgba{1, 0, 0, 0.25},
          .enable_alpha_sliders = true,
      },
      &builder, nullptr, nullptr, nullptr, meshcat);

---

examples/conveyor_belt/conveyor_belt.cc line 54 at r8 (raw file):

  simulator.set_target_realtime_rate(1.0);
  simulator.Initialize();
  visualizer.StartRecording();

Changing visualization requires a change here. visualizer is no longer defined. But it's easily fixed. (Here and on visualizer.PublishRecording();.

Suggestion:

 meshcat->StartRecording();

---

examples/conveyor_belt/conveyor_belt.cc line 71 at r8 (raw file):

  drake::examples::conveyor_belt::do_main_continous_plant();
  return 0;
}

nit: Missing new line.

[juanmed]

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

This badly needs a remerge with master:

1. The build infrastructure has moved on (new version of bazel).
   • Without doing this, it's incredibly cumbersome to build this locally if I have otherwise moved my local drake branch to master and updated my build infrastructure.
2. There are now some heavy conflicts between some of this PR and master: the changes to shape_specification (such that it depends on downstream elements like bvh/obb). Add Aabb-{plane|half space} intersection and missing bindings #24001 introduced conflicting implementations with what has been added in this PR. That implementation explicitly depends on shape specification (by naming HalfSpace).

I have a couple of ideas for testing my assertions about the potential defects in implementation but, in order to verify them, I need to be able to build. So, I'm going to merge master and reconcile conflicts without changing the spirit of the implementation despite the fact that I believe we'll be removing most (if not all) of the geometry changes subsequently.

Hi! It is great to get some attention here. I will slowly reply to your comments.

• Yeah, its been a while, I stopped asking for reviews after several months and is only natural my branch is out of sync
• Ok, if there is an explict dependency on HalfSpace, that will not work with my current implementation.
• If you are able to merge master, sounds good to me. I have no idea what new changes resulted in HalfSpace now being required, nor the bandwidth to get into the weeds. So I guess you are better suited for that. If I find the time to try it, I will let you know.

[SeanCurtis-TRI]

@SeanCurtis-TRI made 1 comment and resolved 1 discussion.
Reviewable status: 34 unresolved discussions, LGTM missing from assignees SeanCurtis-TRI(platform),joemasterjohn, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on juanmed).

---

a discussion (no related file):

Previously, juanmed (Juan) wrote…

Hi! It is great to get some attention here. I will slowly reply to your comments.

• Yeah, its been a while, I stopped asking for reviews after several months and is only natural my branch is out of sync
• Ok, if there is an explict dependency on HalfSpace, that will not work with my current implementation.
• If you are able to merge master, sounds good to me. I have no idea what new changes resulted in HalfSpace now being required, nor the bandwidth to get into the weeds. So I guess you are better suited for that. If I find the time to try it, I will let you know.

Ok, I'll attempt to push an update to the current state that reconciles with current master. (I'm assuming I have the ability to push to your branch). After I do so, it'll be important that you pull from your remote repo before doing any further editing.

On the same vein, some of the comments I've added below I've written some code confirming ideas. If you like, I can push some additional commits as well.

[SeanCurtis-TRI]

I've submitted a few new commits:

1. Merged master (so, we're up to date) -- this includes an additional commit to reconcile the state of the PR with master.
2. Updated the conveyor belt demo to use preferred visualization techniques.
3. Changed how relative velocity is computed. This is the big one.
   • It fixes the defect listed below about creating a relative surface velocity.
   • It no longer assumes the geometry pose is the identity.
   • It also simplifies the computation by using the contact data for normal.

The changes in (3) make all/most of the changes to geometry/ unnecessary. So, we'll want to clean those out. I didn't do so, because I wanted you to see it first. We can also make changes to conveyor_belt.sdf so you can see the difference before and after this third change has been added.

At this point, I'd suggest the next steps:

1. Take a look at the next steps and see if you're persuaded.
2. Assuming you agree, let's back out the now unnecessary changes in geometry.
3. Once clean, we need to revisit the tests to make sure we've truly got the coverage we need.

@SeanCurtis-TRI made 1 comment.
Reviewable status: 34 unresolved discussions, LGTM missing from assignees SeanCurtis-TRI(platform),joemasterjohn, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on juanmed).

[SeanCurtis-TRI]

Somehow the unit test is failing; I'll take a look....

[SeanCurtis-TRI]

I believe I've resolved the unit test problems. Make sure you pull from your forked repo before adding new commits.

[SeanCurtis-TRI]

@drake-jenkins-bot linux-noble-clang-bazel-experimental-everything-release please

[SeanCurtis-TRI]

@juanmed It's been a week, and I thought I'd ping you again. Let me know what you think the likelihood of you coming back to this is. Thanks.

@SeanCurtis-TRI made 1 comment.
Reviewable status: 34 unresolved discussions, LGTM missing from assignees SeanCurtis-TRI(platform),joemasterjohn, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on juanmed).

[SeanCurtis-TRI]

@SeanCurtis-TRI resolved 1 discussion.
Reviewable status: 33 unresolved discussions, LGTM missing from assignees SeanCurtis-TRI(platform),joemasterjohn, commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on juanmed).