@@ -9105,7 +9105,7 @@ Extraction
 
 .. _sensor-contact-data:
 
-:at:`data`: :at-val:`[found, force, torque, dist, pos, normal, tangent], "found"`
+:at:`data`: :at-val:`[found, force, torque, dist, pos, normal, tangent, linvel, angvel], "found"`
    Specification of which data field(s) to report from the selected contacts.
 
    - :at-val:`found` **real(1)**: This field serves two purposes. First, it indicates whether a contact was found in
@@ -9119,6 +9119,11 @@ Extraction
    - :at-val:`normal`: **real(3)**: The contact normal direction, in the global frame.
    - :at-val:`tangent`: **real(3)**: The first tangent direction, in the global frame.
      In order to complete the full 3x3 contact frame, use tangent2 = cross(normal, tangent).
+   - :at-val:`linvel`: **real(3)**: The relative linear velocity of the second contacting body with respect to the first
+     at the contact position, expressed in the contact frame. Positive normal velocity (first component) corresponds to
+     separation. Cannot be combined with ``reduce="netforce"``.
+   - :at-val:`angvel`: **real(3)**: The relative angular velocity of the second contacting body with respect to the
+     first, expressed in the contact frame. Cannot be combined with ``reduce="netforce"``.
 
    Importantly, the :at:`data` attribute can contain **multiple sequential data types**, as long as the relative
    order---as listed above---is maintained. For example, :at:`data` = :at-val:`"found force dist"` will return 5 numbers
@@ -9152,8 +9157,9 @@ Extraction
    Computation stage
       The computation stage (``sensor_needstage``) of a contact sensor depends on the requested :at:`data` fields and
       :ref:`reduce<sensor-contact-reduce>` criterion. If :at:`data` includes ``force`` or ``torque``, or if
-      :at:`reduce` is ``maxforce`` or ``netforce``, the sensor is evaluated at ``mjSTAGE_ACC``; otherwise it is
-      evaluated at ``mjSTAGE_POS``.
+      :at:`reduce` is ``maxforce`` or ``netforce``, the sensor is evaluated at ``mjSTAGE_ACC``; otherwise, if
+      :at:`data` includes ``linvel`` or ``angvel``, it is evaluated at ``mjSTAGE_VEL``; otherwise it is evaluated at
+      ``mjSTAGE_POS``.
 
 .. _sensor-contact-reduce:
 
 
@@ -34,6 +34,9 @@ Engine
   degrees of freedom. Other flexes keep all their contacts. Since version 3.6.0 the contacts of every flex were reduced,
   and in the midphase for all geoms of a body together, so a cloth dropped on a plane or on many static geoms in the
   world body was held by at most 50 contacts and fell through.
+- The :ref:`contact<sensor-contact>` sensor now supports ``linvel`` and ``angvel`` :ref:`data<sensor-contact-data>`
+  fields for reporting relative contact velocities, and its computation stage (:ref:`mjtStage`) is now inferred from the
+  requested :ref:`data<sensor-contact-data>` and :ref:`reduce<sensor-contact-reduce>` attributes.
 - Fixed floating-point cancellation in capsule-capsule collision detection for parallel and near-parallel axes. Thanks
   to :github:user:`ayonga` for reporting the issue in :doc:`MJWarp <mjwarp/index>` issue
   `1758 <https://github.com/google-deepmind/mujoco_warp/issues/1758>`__.
 
@@ -2747,6 +2747,8 @@ typedef enum mjtConDataField {    // data fields returned by contact sensors
   mjCONDATA_POS,                  // contact position
   mjCONDATA_NORMAL,               // contact frame normal
   mjCONDATA_TANGENT,              // contact frame first tangent
+  mjCONDATA_LINVEL,               // contact linear velocity
+  mjCONDATA_ANGVEL,               // contact angular velocity
 
   mjNCONDATA                      // number of contact sensor data fields
 } mjtConDataField;
 
@@ -432,6 +432,8 @@ typedef enum mjtConDataField {    // data fields returned by contact sensors
   mjCONDATA_POS,                  // contact position
   mjCONDATA_NORMAL,               // contact frame normal
   mjCONDATA_TANGENT,              // contact frame first tangent
+  mjCONDATA_LINVEL,               // contact linear velocity
+  mjCONDATA_ANGVEL,               // contact angular velocity
 
   mjNCONDATA                      // number of contact sensor data fields
 } mjtConDataField;
 
@@ -94,9 +94,9 @@ def _reduce(reduction, mask):
   # number of data elements per slot
   def nslotdata(dataspec):
     size = 0
-    # found, force, torque, dist, pos, normal, tangent
+    # found, force, torque, dist, pos, normal, tangent, linvel, angvel
     # TODO(taylorhowell): get sizes from mjCONDATA_SIZE
-    for i, size_i in enumerate([1, 3, 3, 1, 3, 3, 3]):
+    for i, size_i in enumerate([1, 3, 3, 1, 3, 3, 3, 3, 3]):
       if dataspec & (1 << i):
         size += size_i
     return size
@@ -214,6 +214,29 @@ def nslotdata(dataspec):
     if dataspec & (1 << 6):  # tangent
       slot.append(flip[:, 2, None] * d._impl.contact.frame[cid, 1])  # pyrefly: ignore[missing-attribute]
 
+    if dataspec & ((1 << 7) | (1 << 8)):  # linvel or angvel
+      con_pos = d._impl.contact.pos[cid]  # pyrefly: ignore[missing-attribute]
+      con_frame = d._impl.contact.frame[cid]  # pyrefly: ignore[missing-attribute]
+      con_geom = d._impl.contact.geom[cid]  # pyrefly: ignore[missing-attribute]
+      bid0 = jp.array(m.geom_bodyid)[con_geom[:, 0]]
+      bid1 = jp.array(m.geom_bodyid)[con_geom[:, 1]]
+      cvel0 = d.cvel[bid0]
+      cvel1 = d.cvel[bid1]
+
+      if dataspec & (1 << 7):  # linvel
+        root0 = jp.array(m.body_rootid)[bid0]
+        root1 = jp.array(m.body_rootid)[bid1]
+        v0 = cvel0[:, 3:] - jp.cross(con_pos - d.subtree_com[root0], cvel0[:, :3])
+        v1 = cvel1[:, 3:] - jp.cross(con_pos - d.subtree_com[root1], cvel1[:, :3])
+        linvel = jax.vmap(lambda f, v: f @ v)(con_frame, v1 - v0)
+        slot.append(flip * linvel)
+
+      if dataspec & (1 << 8):  # angvel
+        angvel = jax.vmap(lambda f, w: f @ w)(
+            con_frame, cvel1[:, :3] - cvel0[:, :3]
+        )
+        slot.append(flip * angvel)
+
     found = jp.tile(jp.arange(num), nsensor) < jp.repeat(nfound, num)
     sensors.append((found[:, None] * jp.hstack(slot)).reshape(-1))
     adrs.append(
@@ -586,6 +609,11 @@ def _cvel_offset(otype, oid):
     elif sensor_type == SensorType.SUBTREEANGMOM:
       sensor = d._impl.subtree_angmom[objid]  # pyrefly: ignore[missing-attribute]
       adr = (adr[:, None] + np.arange(3)[None]).reshape(-1)
+    elif sensor_type == SensorType.CONTACT:
+      s, a = _sensor_contact(m, d, idx)
+      sensors.extend(s)
+      adrs.extend(a)
+      continue
     else:
       # TODO(taylorhowell): raise error after adding sensor check to io.py
       continue  # unsupported sensor type
 
@@ -196,6 +196,52 @@ def test_sensor_contact(self, geom):
 
     _assert_eq(dx.sensordata, d.sensordata, 'sensordata')
 
+  def test_sensor_contact_velocity(self):
+    """Tests contact sensor linvel and angvel with nonzero velocities."""
+    contact_sensors = ''
+    for data in ['linvel', 'angvel', 'found linvel angvel', 'dist linvel']:
+      for match in [
+          '',
+          'geom1="plane"',
+          'geom1="sphere1"',
+          'geom1="sphere1" geom2="sphere2"',
+          'geom1="sphere2" geom2="sphere1"',
+      ]:
+        for reduce in ['mindist', 'maxforce']:
+          contact_sensors += (
+              f'<contact {match} num="2" data="{data}" reduce="{reduce}"/>\n'
+          )
+
+    m = mujoco.MjModel.from_xml_string(f"""
+      <mujoco>
+        <worldbody>
+          <geom name="plane" type="plane" size="10 10 .001"/>
+          <body pos="0 0 .09">
+            <freejoint/>
+            <geom name="sphere1" type="sphere" size=".1"/>
+          </body>
+          <body pos=".02 0 .27">
+            <freejoint/>
+            <geom name="sphere2" type="sphere" size=".1"/>
+          </body>
+        </worldbody>
+        <sensor>
+          {contact_sensors}
+        </sensor>
+      </mujoco>
+    """)
+    d = mujoco.MjData(m)
+    d.qvel = [.3, -.2, .1, .5, -.4, .7, -.1, .4, -.3, -.6, .2, .8]
+
+    mx = mjx.put_model(m)
+    dx = mjx.put_data(m, d)
+
+    mujoco.mj_forward(m, d)
+    dx = jax.jit(mjx.forward)(mx, dx)
+
+    self.assertTrue(np.any(d.sensordata != 0))
+    _assert_eq(dx.sensordata, d.sensordata, 'sensordata')
+
   def test_unsupported_sensor(self):
     """Tests unsupported sensor raises error."""
     m = mujoco.MjModel.from_xml_string("""
 
@@ -360,6 +360,11 @@ def put_model(mjm: mujoco.MjModel, batch_sizes: dict[str, int] | None = None) ->
   if (mjm.sensor_plugin != -1).any():
     raise NotImplementedError("Sensor plugins not supported.")
 
+  # TODO(team): contact sensor linvel and angvel (mjSTAGE_VEL)
+  contact_dataspec = mjm.sensor_intprm[mjm.sensor_type == mujoco.mjtSensor.mjSENS_CONTACT, 0]
+  if (contact_dataspec & ((1 << 7) | (1 << 8))).any():
+    raise NotImplementedError("Contact sensor linvel and angvel not supported.")
+
   if mjm.nflex > 0:
     for fi in range(mjm.nflex):
       if abs(mjm.flex_interp[fi]) == 2:
 
@@ -448,7 +448,9 @@
              ('mjCONDATA_POS', 4),
              ('mjCONDATA_NORMAL', 5),
              ('mjCONDATA_TANGENT', 6),
-             ('mjNCONDATA', 7),
+             ('mjCONDATA_LINVEL', 7),
+             ('mjCONDATA_ANGVEL', 8),
+             ('mjNCONDATA', 9),
          ]),
      )),
     ('mjtRayDataField',
 
@@ -24,6 +24,7 @@
 #include "engine/engine_callback.h"
 #include "engine/engine_collision_convex.h"
 #include "engine/engine_collision_sdf.h"
+#include "engine/engine_core_constraint.h"
 #include "engine/engine_core_smooth.h"
 #include "engine/engine_core_util.h"
 #include "engine/engine_crossplatform.h"
@@ -392,8 +393,8 @@ static int matchContact(const mjModel* m, const mjData* d, int conid,
 
 // fill in output data for contact sensor for all fields
 //   if flg_flip > 0, normal/tangent rotate 180 about frame[2]
-//   force/torque flip-z s.t. force is equal-and-opposite in new contact frame
-static void copySensorData(const mjModel* m, const mjData* d,
+//   force/torque/linvel/angvel flip-z s.t. values are equal-and-opposite in new contact frame
+static void copySensorData(const mjModel* m, mjData* d,
                            mjtNum* data[mjNCONDATA], int id, int flg_flip, int nfound) {
   // found flag
   if (data[mjCONDATA_FOUND]) *data[mjCONDATA_FOUND] = nfound;
@@ -433,6 +434,65 @@ static void copySensorData(const mjModel* m, const mjData* d,
     mju_copy3(data[mjCONDATA_TANGENT], d->contact[id].frame+3);
     if (flg_flip) mju_scl3(data[mjCONDATA_TANGENT], data[mjCONDATA_TANGENT], -1);
   }
+
+  // contact linear and angular velocity
+  if (data[mjCONDATA_LINVEL] || data[mjCONDATA_ANGVEL]) {
+    const mjContact* con = d->contact + id;
+    mjtNum linvel[3] = {0}, angvel[3] = {0};
+
+    // single body on each side: compute from cvel
+    if ((con->geom[0] >= 0 || (con->vert[0] >= 0 && m->flex_interp[con->flex[0]] == 0)) &&
+        (con->geom[1] >= 0 || (con->vert[1] >= 0 && m->flex_interp[con->flex[1]] == 0))) {
+      int bid[2];
+      mjtNum vel[2][6];
+      for (int s=0; s < 2; s++) {
+        bid[s] = (con->geom[s] >= 0) ?
+                 m->geom_bodyid[con->geom[s]] :
+                 m->flex_vertbodyid[m->flex_vertadr[con->flex[s]] + con->vert[s]];
+        mju_transformSpatial(vel[s], d->cvel + 6*bid[s], 0, con->pos,
+                             d->subtree_com + 3*m->body_rootid[bid[s]], NULL);
+      }
+      mju_sub3(angvel, vel[1], vel[0]);
+      mju_sub3(linvel, vel[1]+3, vel[0]+3);
+    }
+
+    // general flex case: compute via contact Jacobian
+    else if (m->nv) {
+      int nv = m->nv, issparse = mj_isSparse(m);
+      int need_ang = (data[mjCONDATA_ANGVEL] != NULL);
+      mj_markStack(d);
+      mjtNum* jacdifp = mjSTACKALLOC(d, 3*nv, mjtNum);
+      mjtNum* jacdifr = need_ang ? mjSTACKALLOC(d, 3*nv, mjtNum) : NULL;
+      int* chain = issparse ? mjSTACKALLOC(d, nv, int) : NULL;
+      int NV = mj_contactJacobian(m, d, con, need_ang ? 6 : 3, jacdifp, jacdifr,
+                                  NULL, NULL, NULL, NULL, chain);
+      if (NV) {
+        if (issparse) {
+          for (int j=0; j < NV; j++) {
+            mjtNum qv = d->qvel[chain[j]];
+            for (int k=0; k < 3; k++) {
+              linvel[k] += jacdifp[k*NV+j] * qv;
+              if (need_ang) angvel[k] += jacdifr[k*NV+j] * qv;
+            }
+          }
+        } else {
+          mju_mulMatVec(linvel, jacdifp, d->qvel, 3, nv);
+          if (need_ang) mju_mulMatVec(angvel, jacdifr, d->qvel, 3, nv);
+        }
+      }
+      mj_freeStack(d);
+    }
+
+    // rotate to contact frame and apply flip
+    if (data[mjCONDATA_LINVEL]) {
+      mju_mulMatVec3(data[mjCONDATA_LINVEL], con->frame, linvel);
+      if (flg_flip) data[mjCONDATA_LINVEL][2] *= -1;
+    }
+    if (data[mjCONDATA_ANGVEL]) {
+      mju_mulMatVec3(data[mjCONDATA_ANGVEL], con->frame, angvel);
+      if (flg_flip) data[mjCONDATA_ANGVEL][2] *= -1;
+    }
+  }
 }
 
 
@@ -1116,6 +1176,10 @@ static void mj_computeSensorVel(const mjModel* m, mjData* d, int i, mjtNum* sens
     mju_copy3(sensordata, d->subtree_angmom+3*objid);
     break;
 
+  case mjSENS_CONTACT:                                // contact
+    compute_contact_sensor(m, d, i, sensordata);
+    break;
+
   case mjSENS_TACTILE:                                // tactile
     {
       mj_markStack(d);
 
@@ -119,7 +119,9 @@ const int mjCONDATA_SIZE[mjNCONDATA] = {
   1,  // mjCONDATA_DIST
   3,  // mjCONDATA_POS
   3,  // mjCONDATA_NORMAL
-  3   // mjCONDATA_TANGENT
+  3,  // mjCONDATA_TANGENT
+  3,  // mjCONDATA_LINVEL
+  3   // mjCONDATA_ANGVEL
 };