diff --git a/benchmarks/cloth/scene.xml b/benchmarks/cloth/scene.xml
index cce2da213..9905597f9 100644
--- a/benchmarks/cloth/scene.xml
+++ b/benchmarks/cloth/scene.xml
@@ -15,7 +15,7 @@
 
 <mujoco model="Poncho">
   <include file="mannequin.xml"/>
-  <option solver="CG" tolerance="1e-6" jacobian="sparse"/>
+  <option integrator="Euler" solver="CG" viscosity="0.05" tolerance="1e-6"/>
 
   <visual>
     <map force="0.1" zfar="30"/>
@@ -31,10 +31,11 @@
   <worldbody>
     <light directional="false" diffuse=".2 .2 .2" specular="0 0 0" pos="0 0 5" dir="0 0 -1"/>
 
-    <flexcomp name="towel" type="grid" count="30 30 1" spacing="0.05 0.05 0.05"
-              radius="0.02" dim="2" rgba="1 0.5 0.5 1" pos="0 0 2" mass=".1">
-      <edge equality="true"/>
-      <elasticity young="3e3" poisson="0" thickness="1e-2" damping="1e-3" elastic2d="bend"/>
+    <flexcomp name="towel" type="grid" count="30 30 1" spacing="0.06 0.06 0.06"
+              radius="0.01" dim="2" rgba="1 0.5 0.5 1" pos="0 0 2" mass="1">
+      <edge equality="true" damping="0.1"/>
+      <elasticity young="3e5" poisson="0" thickness="8e-3" elastic2d="bend"/>
+      <contact selfcollide="narrow"/>
     </flexcomp>
 
   </worldbody>
diff --git a/mujoco_warp/_src/collision_driver_test.py b/mujoco_warp/_src/collision_driver_test.py
index 6c0397372..225eccbbd 100644
--- a/mujoco_warp/_src/collision_driver_test.py
+++ b/mujoco_warp/_src/collision_driver_test.py
@@ -30,6 +30,7 @@
 from mujoco_warp._src.collision_driver import MJ_COLLISION_TABLE
 from mujoco_warp._src.collision_primitive import plane_convex
 from mujoco_warp._src.math import upper_trid_index
+from mujoco_warp._src.util_pkg import check_version
 from mujoco_warp.test_data.collision_sdf.utils import register_sdf_plugins
 
 _TOLERANCE = 5e-5
@@ -572,7 +573,10 @@ def test_collision(self, fixture):
           break
       np.testing.assert_equal(result, True, f"Contact {i} not found in Gjk results")
 
-    self.assertEqual(d.nacon.numpy()[0], mjd.ncon)
+      if check_version("mujoco>=3.5.1.dev888028695"):
+        self.assertGreaterEqual(d.nacon.numpy()[0], mjd.ncon)
+      else:
+        self.assertEqual(d.nacon.numpy()[0], mjd.ncon)
 
   _HFIELD_FIXTURES = {
     "hfield_box": """
diff --git a/mujoco_warp/_src/collision_flex.py b/mujoco_warp/_src/collision_flex.py
index 13d8b0e0d..4ba002142 100644
--- a/mujoco_warp/_src/collision_flex.py
+++ b/mujoco_warp/_src/collision_flex.py
@@ -19,7 +19,9 @@
 from mujoco_warp._src.math import make_frame
 from mujoco_warp._src.types import MJ_MAXVAL
 from mujoco_warp._src.types import MJ_MINMU
+from mujoco_warp._src.types import ContactType
 from mujoco_warp._src.types import Data
+from mujoco_warp._src.types import FlexSelfCollideType
 from mujoco_warp._src.types import GeomType
 from mujoco_warp._src.types import Model
 from mujoco_warp._src.types import vec5
@@ -41,9 +43,9 @@ def _write_flex_contact(
   friction: vec5,
   solref: wp.vec2,
   solimp: vec5,
-  geom: int,
-  flexid: int,
-  vertid: int,
+  geom_ids: wp.vec2i,
+  flex_ids: wp.vec2i,
+  vert_ids: wp.vec2i,
   worldid: int,
   # Data out:
   contact_dist_out: wp.array(dtype=float),
@@ -60,7 +62,6 @@ def _write_flex_contact(
   contact_vert_out: wp.array(dtype=wp.vec2i),
   contact_worldid_out: wp.array(dtype=int),
   contact_type_out: wp.array(dtype=int),
-  contact_geomcollisionid_out: wp.array(dtype=int),
   nacon_out: wp.array(dtype=int),
 ):
   if dist >= margin or dist >= MJ_MAXVAL:
@@ -79,12 +80,11 @@ def _write_flex_contact(
   contact_solreffriction_out[id_] = wp.vec2(0.0, 0.0)
   contact_solimp_out[id_] = solimp
   contact_dim_out[id_] = condim
-  contact_geom_out[id_] = wp.vec2i(geom, -1)
-  contact_flex_out[id_] = wp.vec2i(-1, flexid)
-  contact_vert_out[id_] = wp.vec2i(-1, vertid)
+  contact_geom_out[id_] = geom_ids
+  contact_flex_out[id_] = flex_ids
+  contact_vert_out[id_] = vert_ids
   contact_worldid_out[id_] = worldid
-  contact_type_out[id_] = 1
-  contact_geomcollisionid_out[id_] = 0
+  contact_type_out[id_] = int(ContactType.CONSTRAINT)
 
 
 @wp.func
@@ -105,9 +105,9 @@ def _collide_geom_triangle(
   friction: vec5,
   solref: wp.vec2,
   solimp: vec5,
-  geomid: int,
-  flexid: int,
-  vertex_id: int,
+  geom_ids: wp.vec2i,
+  flex_ids: wp.vec2i,
+  vert_ids: wp.vec2i,
   worldid: int,
   # Data out:
   contact_dist_out: wp.array(dtype=float),
@@ -124,7 +124,6 @@ def _collide_geom_triangle(
   contact_vert_out: wp.array(dtype=wp.vec2i),
   contact_worldid_out: wp.array(dtype=int),
   contact_type_out: wp.array(dtype=int),
-  contact_geomcollisionid_out: wp.array(dtype=int),
   nacon_out: wp.array(dtype=int),
 ):
   if gtype == int(GeomType.SPHERE):
@@ -141,9 +140,9 @@ def _collide_geom_triangle(
         friction,
         solref,
         solimp,
-        geomid,
-        flexid,
-        vertex_id,
+        geom_ids,
+        flex_ids,
+        vert_ids,
         worldid,
         contact_dist_out,
         contact_pos_out,
@@ -159,12 +158,10 @@ def _collide_geom_triangle(
         contact_vert_out,
         contact_worldid_out,
         contact_type_out,
-        contact_geomcollisionid_out,
         nacon_out,
       )
     return
 
-  # Capsule, box, cylinder all return up to 2 contacts - compute then share writing code
   dists = wp.vec2(collision_primitive_core.MJ_MAXVAL, collision_primitive_core.MJ_MAXVAL)
   poss = collision_primitive_core.mat23f(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
   nrms = collision_primitive_core.mat23f(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
@@ -186,7 +183,6 @@ def _collide_geom_triangle(
       pos, cyl_axis, cyl_radius, cyl_half_height, t1, t2, t3, tri_radius
     )
 
-  # Write up to 2 contacts (shared code for capsule/box/cylinder)
   if dists[0] < margin:
     p1 = wp.vec3(poss[0, 0], poss[0, 1], poss[0, 2])
     n1 = wp.vec3(nrms[0, 0], nrms[0, 1], nrms[0, 2])
@@ -200,9 +196,9 @@ def _collide_geom_triangle(
       friction,
       solref,
       solimp,
-      geomid,
-      flexid,
-      vertex_id,
+      geom_ids,
+      flex_ids,
+      vert_ids,
       worldid,
       contact_dist_out,
       contact_pos_out,
@@ -218,7 +214,6 @@ def _collide_geom_triangle(
       contact_vert_out,
       contact_worldid_out,
       contact_type_out,
-      contact_geomcollisionid_out,
       nacon_out,
     )
   if dists[1] < margin:
@@ -234,9 +229,9 @@ def _collide_geom_triangle(
       friction,
       solref,
       solimp,
-      geomid,
-      flexid,
-      vertex_id,
+      geom_ids,
+      flex_ids,
+      vert_ids,
       worldid,
       contact_dist_out,
       contact_pos_out,
@@ -252,7 +247,6 @@ def _collide_geom_triangle(
       contact_vert_out,
       contact_worldid_out,
       contact_type_out,
-      contact_geomcollisionid_out,
       nacon_out,
     )
 
@@ -295,7 +289,6 @@ def _flex_plane_narrowphase(
   contact_vert_out: wp.array(dtype=wp.vec2i),
   contact_worldid_out: wp.array(dtype=int),
   contact_type_out: wp.array(dtype=int),
-  contact_geomcollisionid_out: wp.array(dtype=int),
   nacon_out: wp.array(dtype=int),
 ):
   worldid, vertid = wp.tid()
@@ -354,9 +347,9 @@ def _flex_plane_narrowphase(
         friction,
         solref,
         solimp,
-        geomid,
-        flexid,
-        local_vertid,
+        wp.vec2i(geomid, -1),
+        wp.vec2i(-1, flexid),
+        wp.vec2i(-1, local_vertid),
         worldid,
         contact_dist_out,
         contact_pos_out,
@@ -372,7 +365,6 @@ def _flex_plane_narrowphase(
         contact_vert_out,
         contact_worldid_out,
         contact_type_out,
-        contact_geomcollisionid_out,
         nacon_out,
       )
 
@@ -393,12 +385,13 @@ def _flex_narrowphase_dim2(
   geom_margin: wp.array2d(dtype=float),
   flex_contype: wp.array(dtype=int),
   flex_conaffinity: wp.array(dtype=int),
+  flex_condim: wp.array(dtype=int),
+  flex_friction: wp.array(dtype=wp.vec3),
   flex_margin: wp.array(dtype=float),
   flex_dim: wp.array(dtype=int),
   flex_vertadr: wp.array(dtype=int),
   flex_elemadr: wp.array(dtype=int),
   flex_elemnum: wp.array(dtype=int),
-  flex_elemdataadr: wp.array(dtype=int),
   flex_elem: wp.array(dtype=int),
   flex_radius: wp.array(dtype=float),
   # Data in:
@@ -422,7 +415,6 @@ def _flex_narrowphase_dim2(
   contact_vert_out: wp.array(dtype=wp.vec2i),
   contact_worldid_out: wp.array(dtype=int),
   contact_type_out: wp.array(dtype=int),
-  contact_geomcollisionid_out: wp.array(dtype=int),
   nacon_out: wp.array(dtype=int),
 ):
   worldid, elemid = wp.tid()
@@ -444,7 +436,7 @@ def _flex_narrowphase_dim2(
   tri_radius = flex_radius[flexid]
   tri_margin = flex_margin[flexid]
 
-  elem_data_idx = flex_elemdataadr[flexid] + (elemid - flex_elemadr[flexid]) * 3
+  elem_data_idx = elemid * 3
   v0_local = flex_elem[elem_data_idx]
   v1_local = flex_elem[elem_data_idx + 1]
   v2_local = flex_elem[elem_data_idx + 2]
@@ -478,14 +470,18 @@ def _flex_narrowphase_dim2(
     geom_rot = geom_xmat_in[worldid, geomid]
     geom_size_val = geom_size[worldid % geom_size.shape[0], geomid]
 
-    condim = geom_condim[geomid]
+    condim = wp.max(geom_condim[geomid], flex_condim[flexid])
     gf = geom_friction[worldid % geom_friction.shape[0], geomid]
+    ff = flex_friction[flexid]
+    fric0 = wp.max(gf[0], ff[0])
+    fric1 = wp.max(gf[1], ff[1])
+    fric2 = wp.max(gf[2], ff[2])
     friction = vec5(
-      wp.max(MJ_MINMU, gf[0]),
-      wp.max(MJ_MINMU, gf[0]),
-      wp.max(MJ_MINMU, gf[1]),
-      wp.max(MJ_MINMU, gf[2]),
-      wp.max(MJ_MINMU, gf[2]),
+      wp.max(MJ_MINMU, fric0),
+      wp.max(MJ_MINMU, fric0),
+      wp.max(MJ_MINMU, fric1),
+      wp.max(MJ_MINMU, fric2),
+      wp.max(MJ_MINMU, fric2),
     )
     solref = geom_solref[worldid % geom_solref.shape[0], geomid]
     solimp = geom_solimp[worldid % geom_solimp.shape[0], geomid]
@@ -505,9 +501,9 @@ def _flex_narrowphase_dim2(
       friction,
       solref,
       solimp,
-      geomid,
-      flexid,
-      v0_local,
+      wp.vec2i(geomid, -1),
+      wp.vec2i(-1, flexid),
+      wp.vec2i(-1, v0_local),
       worldid,
       contact_dist_out,
       contact_pos_out,
@@ -523,7 +519,6 @@ def _flex_narrowphase_dim2(
       contact_vert_out,
       contact_worldid_out,
       contact_type_out,
-      contact_geomcollisionid_out,
       nacon_out,
     )
 
@@ -544,6 +539,8 @@ def _flex_narrowphase_dim3(
   geom_margin: wp.array2d(dtype=float),
   flex_contype: wp.array(dtype=int),
   flex_conaffinity: wp.array(dtype=int),
+  flex_condim: wp.array(dtype=int),
+  flex_friction: wp.array(dtype=wp.vec3),
   flex_margin: wp.array(dtype=float),
   flex_dim: wp.array(dtype=int),
   flex_vertadr: wp.array(dtype=int),
@@ -572,7 +569,6 @@ def _flex_narrowphase_dim3(
   contact_vert_out: wp.array(dtype=wp.vec2i),
   contact_worldid_out: wp.array(dtype=int),
   contact_type_out: wp.array(dtype=int),
-  contact_geomcollisionid_out: wp.array(dtype=int),
   nacon_out: wp.array(dtype=int),
 ):
   worldid, shellid = wp.tid()
@@ -632,14 +628,18 @@ def _flex_narrowphase_dim3(
     geom_rot = geom_xmat_in[worldid, geomid]
     geom_size_val = geom_size[worldid % geom_size.shape[0], geomid]
 
-    condim = geom_condim[geomid]
+    condim = wp.max(geom_condim[geomid], flex_condim[flexid])
     gf = geom_friction[worldid % geom_friction.shape[0], geomid]
+    ff = flex_friction[flexid]
+    fric0 = wp.max(gf[0], ff[0])
+    fric1 = wp.max(gf[1], ff[1])
+    fric2 = wp.max(gf[2], ff[2])
     friction = vec5(
-      wp.max(MJ_MINMU, gf[0]),
-      wp.max(MJ_MINMU, gf[0]),
-      wp.max(MJ_MINMU, gf[1]),
-      wp.max(MJ_MINMU, gf[2]),
-      wp.max(MJ_MINMU, gf[2]),
+      wp.max(MJ_MINMU, fric0),
+      wp.max(MJ_MINMU, fric0),
+      wp.max(MJ_MINMU, fric1),
+      wp.max(MJ_MINMU, fric2),
+      wp.max(MJ_MINMU, fric2),
     )
     solref = geom_solref[worldid % geom_solref.shape[0], geomid]
     solimp = geom_solimp[worldid % geom_solimp.shape[0], geomid]
@@ -659,9 +659,212 @@ def _flex_narrowphase_dim3(
       friction,
       solref,
       solimp,
-      geomid,
-      flexid,
-      v0_local,
+      wp.vec2i(geomid, -1),
+      wp.vec2i(-1, flexid),
+      wp.vec2i(-1, v0_local),
+      worldid,
+      contact_dist_out,
+      contact_pos_out,
+      contact_frame_out,
+      contact_includemargin_out,
+      contact_friction_out,
+      contact_solref_out,
+      contact_solreffriction_out,
+      contact_solimp_out,
+      contact_dim_out,
+      contact_geom_out,
+      contact_flex_out,
+      contact_vert_out,
+      contact_worldid_out,
+      contact_type_out,
+      nacon_out,
+    )
+
+
+@wp.func
+def _share_vertex_inline(
+  # Model:
+  flex_elem: wp.array(dtype=int),
+  # In:
+  e1_global: int,
+  e2_global: int,
+  dim: int,
+) -> int:
+  for i in range(dim + 1):
+    v1 = flex_elem[e1_global * (dim + 1) + i]
+    for j in range(dim + 1):
+      v2 = flex_elem[e2_global * (dim + 1) + j]
+      if v1 == v2:
+        return 1
+  return 0
+
+
+@wp.func
+def _aabb_overlap_3d(
+  # In:
+  aabb: wp.array3d(dtype=float),
+  worldid: int,
+  e1: int,
+  e2: int,
+) -> int:
+  for axis in range(3):
+    c1 = aabb[worldid, e1, axis]
+    r1 = aabb[worldid, e1, axis + 3]
+    c2 = aabb[worldid, e2, axis]
+    r2 = aabb[worldid, e2, axis + 3]
+    if wp.abs(c1 - c2) > r1 + r2:
+      return 0
+  return 1
+
+
+@wp.kernel
+def _flex_self_narrow(
+  # Model:
+  nflex: int,
+  flex_condim: wp.array(dtype=int),
+  flex_solref: wp.array(dtype=wp.vec2),
+  flex_solimp: wp.array(dtype=vec5),
+  flex_friction: wp.array(dtype=wp.vec3),
+  flex_selfcollide: wp.array(dtype=int),
+  flex_dim: wp.array(dtype=int),
+  flex_vertadr: wp.array(dtype=int),
+  flex_elemadr: wp.array(dtype=int),
+  flex_elemnum: wp.array(dtype=int),
+  flex_elem: wp.array(dtype=int),
+  flex_radius: wp.array(dtype=float),
+  # Data in:
+  flexvert_xpos_in: wp.array2d(dtype=wp.vec3),
+  flexelem_aabb_in: wp.array3d(dtype=float),
+  naconmax_in: int,
+  # Data out:
+  contact_dist_out: wp.array(dtype=float),
+  contact_pos_out: wp.array(dtype=wp.vec3),
+  contact_frame_out: wp.array(dtype=wp.mat33),
+  contact_includemargin_out: wp.array(dtype=float),
+  contact_friction_out: wp.array(dtype=vec5),
+  contact_solref_out: wp.array(dtype=wp.vec2),
+  contact_solreffriction_out: wp.array(dtype=wp.vec2),
+  contact_solimp_out: wp.array(dtype=vec5),
+  contact_dim_out: wp.array(dtype=int),
+  contact_geom_out: wp.array(dtype=wp.vec2i),
+  contact_flex_out: wp.array(dtype=wp.vec2i),
+  contact_vert_out: wp.array(dtype=wp.vec2i),
+  contact_worldid_out: wp.array(dtype=int),
+  contact_type_out: wp.array(dtype=int),
+  nacon_out: wp.array(dtype=int),
+):
+  worldid, pairid = wp.tid()
+
+  for f in range(nflex):
+    if flex_selfcollide[f] == int(FlexSelfCollideType.NONE):
+      continue
+    if flex_dim[f] != 2:
+      continue
+
+    nelem = flex_elemnum[f]
+    npairs = nelem * (nelem - 1) // 2
+
+    if pairid >= npairs:
+      continue
+
+    i = int(wp.floor((wp.sqrt(8.0 * float(pairid) + 1.0) - 1.0) / 2.0))
+    while i * (i + 1) // 2 <= pairid:
+      i += 1
+    e2 = i
+    e1 = pairid - e2 * (e2 - 1) // 2
+
+    if e1 >= e2 or e2 >= nelem:
+      continue
+
+    dim = flex_dim[f]
+    elemadr = flex_elemadr[f]
+    e1_global = elemadr + e1
+    e2_global = elemadr + e2
+    vertadr = flex_vertadr[f]
+
+    shared = _share_vertex_inline(flex_elem, e1_global, e2_global, dim)
+    if shared:
+      continue
+
+    aabb_ok = _aabb_overlap_3d(flexelem_aabb_in, worldid, e1_global, e2_global)
+    if not aabb_ok:
+      continue
+
+    v1_0 = flex_elem[e1_global * 3 + 0]
+    v1_1 = flex_elem[e1_global * 3 + 1]
+    v1_2 = flex_elem[e1_global * 3 + 2]
+    v2_0 = flex_elem[e2_global * 3 + 0]
+    v2_1 = flex_elem[e2_global * 3 + 1]
+    v2_2 = flex_elem[e2_global * 3 + 2]
+
+    p1_0 = flexvert_xpos_in[worldid, vertadr + v1_0]
+    p1_1 = flexvert_xpos_in[worldid, vertadr + v1_1]
+    p1_2 = flexvert_xpos_in[worldid, vertadr + v1_2]
+    p2_0 = flexvert_xpos_in[worldid, vertadr + v2_0]
+    p2_1 = flexvert_xpos_in[worldid, vertadr + v2_1]
+    p2_2 = flexvert_xpos_in[worldid, vertadr + v2_2]
+
+    radius = flex_radius[f]
+
+    dist, pos, normal = collision_primitive_core.triangle_triangle(
+      p1_0,
+      p1_1,
+      p1_2,
+      radius,
+      p2_0,
+      p2_1,
+      p2_2,
+      radius,
+    )
+
+    margin = 0.0
+    frame = make_frame(normal)
+
+    condim = flex_condim[f]
+    fric = flex_friction[f]
+    friction = vec5(
+      wp.max(MJ_MINMU, fric[0]),
+      wp.max(MJ_MINMU, fric[0]),
+      wp.max(MJ_MINMU, fric[1]),
+      wp.max(MJ_MINMU, fric[2]),
+      wp.max(MJ_MINMU, fric[2]),
+    )
+    solref = flex_solref[f]
+    solimp = flex_solimp[f]
+
+    best_v1 = v1_0
+    best_d1 = wp.length(pos - p1_0)
+    d_tmp = wp.length(pos - p1_1)
+    if d_tmp < best_d1:
+      best_d1 = d_tmp
+      best_v1 = v1_1
+    d_tmp = wp.length(pos - p1_2)
+    if d_tmp < best_d1:
+      best_v1 = v1_2
+
+    best_v2 = v2_0
+    best_d2 = wp.length(pos - p2_0)
+    d_tmp = wp.length(pos - p2_1)
+    if d_tmp < best_d2:
+      best_d2 = d_tmp
+      best_v2 = v2_1
+    d_tmp = wp.length(pos - p2_2)
+    if d_tmp < best_d2:
+      best_v2 = v2_2
+
+    _write_flex_contact(
+      naconmax_in,
+      dist,
+      pos,
+      frame,
+      margin,
+      condim,
+      friction,
+      solref,
+      solimp,
+      wp.vec2i(-1, -1),
+      wp.vec2i(f, f),
+      wp.vec2i(best_v1, best_v2),
       worldid,
       contact_dist_out,
       contact_pos_out,
@@ -677,7 +880,6 @@ def _flex_narrowphase_dim3(
       contact_vert_out,
       contact_worldid_out,
       contact_type_out,
-      contact_geomcollisionid_out,
       nacon_out,
     )
 
@@ -705,12 +907,13 @@ def flex_narrowphase(m: Model, d: Data):
       m.geom_margin,
       m.flex_contype,
       m.flex_conaffinity,
+      m.flex_condim,
+      m.flex_friction,
       m.flex_margin,
       m.flex_dim,
       m.flex_vertadr,
       m.flex_elemadr,
       m.flex_elemnum,
-      m.flex_elemdataadr,
       m.flex_elem,
       m.flex_radius,
       d.geom_xpos,
@@ -734,7 +937,6 @@ def flex_narrowphase(m: Model, d: Data):
       d.contact.vert,
       d.contact.worldid,
       d.contact.type,
-      d.contact.geomcollisionid,
       d.nacon,
     ],
   )
@@ -756,6 +958,8 @@ def flex_narrowphase(m: Model, d: Data):
       m.geom_margin,
       m.flex_contype,
       m.flex_conaffinity,
+      m.flex_condim,
+      m.flex_friction,
       m.flex_margin,
       m.flex_dim,
       m.flex_vertadr,
@@ -784,7 +988,6 @@ def flex_narrowphase(m: Model, d: Data):
       d.contact.vert,
       d.contact.worldid,
       d.contact.type,
-      d.contact.geomcollisionid,
       d.nacon,
     ],
   )
@@ -828,7 +1031,49 @@ def flex_narrowphase(m: Model, d: Data):
       d.contact.vert,
       d.contact.worldid,
       d.contact.type,
-      d.contact.geomcollisionid,
       d.nacon,
     ],
   )
+
+  # Self-collision: only for dim=2 elements with selfcollide enabled.
+  # flex_self_max_pairs is precomputed in put_model to avoid device→host
+  # copies during CUDA graph capture.
+  if m.flex_self_max_pairs > 0:
+    wp.launch(
+      _flex_self_narrow,
+      dim=(d.nworld, m.flex_self_max_pairs),
+      inputs=[
+        m.nflex,
+        m.flex_condim,
+        m.flex_solref,
+        m.flex_solimp,
+        m.flex_friction,
+        m.flex_selfcollide,
+        m.flex_dim,
+        m.flex_vertadr,
+        m.flex_elemadr,
+        m.flex_elemnum,
+        m.flex_elem,
+        m.flex_radius,
+        d.flexvert_xpos,
+        d.flexelem_aabb,
+        d.naconmax,
+      ],
+      outputs=[
+        d.contact.dist,
+        d.contact.pos,
+        d.contact.frame,
+        d.contact.includemargin,
+        d.contact.friction,
+        d.contact.solref,
+        d.contact.solreffriction,
+        d.contact.solimp,
+        d.contact.dim,
+        d.contact.geom,
+        d.contact.flex,
+        d.contact.vert,
+        d.contact.worldid,
+        d.contact.type,
+        d.nacon,
+      ],
+    )
diff --git a/mujoco_warp/_src/collision_flex_test.py b/mujoco_warp/_src/collision_flex_test.py
index 2f61f8e5f..8434e18f3 100644
--- a/mujoco_warp/_src/collision_flex_test.py
+++ b/mujoco_warp/_src/collision_flex_test.py
@@ -13,6 +13,9 @@
 # limitations under the License.
 """Tests for flex element collision."""
 
+import mujoco
+import numpy as np
+import warp as wp
 from absl.testing import absltest
 
 import mujoco_warp as mjwarp
@@ -67,6 +70,141 @@ def test_sphere_cloth_contact_generated(self):
     # Sphere is just above the cloth, so there should be contacts
     self.assertGreater(nacon, 0, "Expected contacts between sphere and cloth")
 
+  def test_folded_cloth_self_contacts_match_c(self):
+    """Folded cloth self-contacts should match between C and Warp."""
+    xml = """
+    <mujoco>
+      <option solver="CG" tolerance="1e-6" timestep=".002"/>
+      <size memory="10M"/>
+      <worldbody>
+        <light pos="0 0 3" dir="0 0 -1"/>
+        <geom type="plane" size="5 5 .1" pos="0 0 -1"/>
+        <flexcomp type="grid" count="8 4 1" spacing=".05 .05 .1"
+                  pos="-.175 -.075 0.5" radius=".005" name="cloth" dim="2"
+                  mass=".5">
+          <contact condim="3" solref="0.01 1" solimp=".95 .99 .0001"
+                   selfcollide="narrow" conaffinity="1" contype="1"/>
+          <edge damping="0.01"/>
+        </flexcomp>
+      </worldbody>
+    </mujoco>
+    """
+    mjm = mujoco.MjModel.from_xml_string(xml)
+    mjd = mujoco.MjData(mjm)
+    mujoco.mj_resetData(mjm, mjd)
+    mujoco.mj_forward(mjm, mjd)
+
+    # Fold the right half of the cloth over the left half.
+    nvert = mjm.flex_vertnum[0]
+    vertadr = mjm.flex_vertadr[0]
+    verts = mjd.flexvert_xpos[vertadr : vertadr + nvert].copy()
+    mid_x = (verts[:, 0].min() + verts[:, 0].max()) / 2.0
+
+    for i in range(nvert):
+      bodyid = mjm.flex_vertbodyid[i]
+      jntid = mjm.body_jntadr[bodyid]
+      qposadr = mjm.jnt_qposadr[jntid]
+      if verts[i, 0] > mid_x + 1e-6:
+        dx = verts[i, 0] - mid_x
+        new_x = mid_x - dx - 0.05
+        delta_x = new_x - verts[i, 0]
+        mjd.qpos[qposadr + 0] = delta_x
+        mjd.qpos[qposadr + 1] = 0.0
+        mjd.qpos[qposadr + 2] = -0.002
+
+    mujoco.mj_forward(mjm, mjd)
+
+    # Count C self-contacts.
+    c_self = 0
+    for c in range(mjd.ncon):
+      con = mjd.contact[c]
+      if con.flex[0] >= 0 and con.flex[0] == con.flex[1]:
+        c_self += 1
+
+    self.assertGreater(c_self, 0, f"C should produce self-contacts")
+
+    # Run Warp with matched state.
+    m = mjwarp.put_model(mjm)
+    d = mjwarp.put_data(mjm, mjd)
+    mjwarp.kinematics(m, d)
+    mjwarp.collision(m, d)
+    wp.synchronize_device()
+
+    nacon = int(d.nacon.numpy()[0])
+    if nacon > 0:
+      flex_arr = d.contact.flex.numpy()[:nacon]
+      warp_self = int(np.sum((flex_arr[:, 0] >= 0) & (flex_arr[:, 0] == flex_arr[:, 1])))
+    else:
+      warp_self = 0
+
+    self.assertGreater(warp_self, 0, f"Warp should produce self-contacts (nacon={nacon})")
+    self.assertGreaterEqual(warp_self, c_self, f"Warp={warp_self} should be >= C={c_self} self-contacts")
+
+  def test_flat_cloth_self_contacts_match_c(self):
+    """Flat cloth (no fold) self-contacts should match C MuJoCo."""
+    xml = """
+    <mujoco>
+      <option solver="CG" tolerance="1e-6" timestep=".002"/>
+      <size memory="10M"/>
+      <worldbody>
+        <flexcomp type="grid" count="10 10 1" spacing=".05 .05 .05"
+                  radius=".02" name="cloth" dim="2" mass=".1"
+                  pos="0 0 2">
+          <contact selfcollide="narrow" condim="3"/>
+          <edge damping="0.01"/>
+        </flexcomp>
+      </worldbody>
+    </mujoco>
+    """
+    mjm_c = mujoco.MjModel.from_xml_string(xml)
+    mjd_c = mujoco.MjData(mjm_c)
+    mujoco.mj_forward(mjm_c, mjd_c)
+    c_self = 0
+    for c in range(mjd_c.ncon):
+      con = mjd_c.contact[c]
+      if con.flex[0] >= 0 and con.flex[0] == con.flex[1]:
+        c_self += 1
+
+    mjd_clean = mujoco.MjData(mjm_c)
+    mujoco.mj_resetData(mjm_c, mjd_clean)
+    mujoco.mj_kinematics(mjm_c, mjd_clean)
+    mujoco.mj_comPos(mjm_c, mjd_clean)
+    m = mjwarp.put_model(mjm_c)
+    d = mjwarp.put_data(mjm_c, mjd_clean)
+    mjwarp.kinematics(m, d)
+    mjwarp.collision(m, d)
+    wp.synchronize_device()
+
+    nacon = int(d.nacon.numpy()[0])
+    if nacon > 0:
+      flex_arr = d.contact.flex.numpy()[:nacon]
+      warp_self = int(np.sum((flex_arr[:, 0] >= 0) & (flex_arr[:, 0] == flex_arr[:, 1])))
+    else:
+      warp_self = 0
+
+    self.assertGreaterEqual(warp_self, c_self, f"Warp={warp_self} should be >= C={c_self} self-contacts")
+
+  def test_unsupported_selfcollide_raises(self):
+    """Verify unsupported selfcollide modes raise NotImplementedError."""
+    xml = """
+    <mujoco>
+      <option solver="CG" tolerance="1e-6" timestep=".002"/>
+      <size memory="10M"/>
+      <worldbody>
+        <flexcomp type="grid" count="3 3 1" spacing=".1 .1 .1" pos="0 0 0"
+                  radius=".01" name="cloth" dim="2" mass=".5">
+          <contact selfcollide="auto" conaffinity="1" contype="1"/>
+          <edge damping="0.01"/>
+        </flexcomp>
+      </worldbody>
+    </mujoco>
+    """
+    mjm = mujoco.MjModel.from_xml_string(xml)
+    mjd = mujoco.MjData(mjm)
+    mujoco.mj_forward(mjm, mjd)
+    with self.assertRaises(NotImplementedError):
+      mjwarp.put_model(mjm)
+
 
 if __name__ == "__main__":
   absltest.main()
diff --git a/mujoco_warp/_src/collision_primitive_core.py b/mujoco_warp/_src/collision_primitive_core.py
index 8a85d2412..0f34fbdec 100644
--- a/mujoco_warp/_src/collision_primitive_core.py
+++ b/mujoco_warp/_src/collision_primitive_core.py
@@ -1596,6 +1596,176 @@ def sphere_triangle(
   return dist, pos, nrm
 
 
+@wp.func
+def _seg_seg_dist(
+  p1: wp.vec3,
+  d1: wp.vec3,
+  p2: wp.vec3,
+  d2: wp.vec3,
+) -> Tuple[float, wp.vec3, wp.vec3]:
+  r = p1 - p2
+  a = wp.dot(d1, d1)
+  e = wp.dot(d2, d2)
+  f = wp.dot(d2, r)
+
+  eps = float(1e-10)
+  s = float(0.0)
+  t = float(0.0)
+  c1 = p1
+  c2 = p2
+
+  if a <= eps and e <= eps:
+    c1 = p1
+    c2 = p2
+  elif a <= eps:
+    t = wp.clamp(f / e, 0.0, 1.0)
+    c1 = p1
+    c2 = p2 + d2 * t
+  elif e <= eps:
+    b = wp.dot(d1, r)
+    s = wp.clamp(-b / a, 0.0, 1.0)
+    c1 = p1 + d1 * s
+    c2 = p2
+  else:
+    b = wp.dot(d1, r)
+    c = wp.dot(d1, d2)
+    denom = a * e - c * c
+    if denom > eps:
+      s = wp.clamp((c * f - b * e) / denom, 0.0, 1.0)
+    else:
+      s = 0.0
+    t = (c * s + f) / e
+    if t < 0.0:
+      t = 0.0
+      s = wp.clamp(-b / a, 0.0, 1.0)
+    elif t > 1.0:
+      t = 1.0
+      s = wp.clamp((c - b) / a, 0.0, 1.0)
+    c1 = p1 + d1 * s
+    c2 = p2 + d2 * t
+
+  return wp.length(c1 - c2), c1, c2
+
+
+@wp.func
+def triangle_triangle(
+  # In:
+  a1: wp.vec3,
+  a2: wp.vec3,
+  a3: wp.vec3,
+  a_radius: float,
+  b1: wp.vec3,
+  b2: wp.vec3,
+  b3: wp.vec3,
+  b_radius: float,
+) -> Tuple[float, wp.vec3, wp.vec3]:
+  best_d = MJ_MAXVAL
+  best_c1 = wp.vec3(0.0)
+  best_c2 = wp.vec3(0.0)
+
+  ea = wp.vec3(0.0)
+  eb = wp.vec3(0.0)
+  pa = wp.vec3(0.0)
+  pb = wp.vec3(0.0)
+  for ia in range(3):
+    if ia == 0:
+      pa = a1
+      ea = a2 - a1
+    elif ia == 1:
+      pa = a2
+      ea = a3 - a2
+    else:
+      pa = a3
+      ea = a1 - a3
+    for ib in range(3):
+      if ib == 0:
+        pb = b1
+        eb = b2 - b1
+      elif ib == 1:
+        pb = b2
+        eb = b3 - b2
+      else:
+        pb = b3
+        eb = b1 - b3
+      d, c1, c2 = _seg_seg_dist(pa, ea, pb, eb)
+      if d < best_d:
+        best_d = d
+        best_c1 = c1
+        best_c2 = c2
+
+  for iv in range(6):
+    vert = a1
+    t1 = b1
+    t2 = b2
+    t3 = b3
+    if iv == 1:
+      vert = a2
+    elif iv == 2:
+      vert = a3
+    elif iv == 3:
+      vert = b1
+      t1 = a1
+      t2 = a2
+      t3 = a3
+    elif iv == 4:
+      vert = b2
+      t1 = a1
+      t2 = a2
+      t3 = a3
+    elif iv == 5:
+      vert = b3
+      t1 = a1
+      t2 = a2
+      t3 = a3
+
+    e1 = t2 - t1
+    e2 = t3 - t1
+    nrm = wp.cross(e1, e2)
+    nrm_len = wp.length(nrm)
+    if nrm_len > 1e-10:
+      nrm = nrm / nrm_len
+      diff = vert - t1
+      proj = vert - nrm * wp.dot(diff, nrm)
+
+      v0 = t3 - t1
+      v1 = t2 - t1
+      v2 = proj - t1
+      d00 = wp.dot(v0, v0)
+      d01 = wp.dot(v0, v1)
+      d02 = wp.dot(v0, v2)
+      d11 = wp.dot(v1, v1)
+      d12 = wp.dot(v1, v2)
+      inv = 1.0 / (d00 * d11 - d01 * d01 + 1e-30)
+      u = (d11 * d02 - d01 * d12) * inv
+      v = (d00 * d12 - d01 * d02) * inv
+
+      if u >= 0.0 and v >= 0.0 and u + v <= 1.0:
+        d = wp.length(vert - proj)
+        if d < best_d:
+          best_d = d
+          if iv < 3:
+            best_c1 = vert
+            best_c2 = proj
+          else:
+            best_c1 = proj
+            best_c2 = vert
+
+  total_radius = a_radius + b_radius
+  dist = best_d - total_radius
+
+  diff = best_c2 - best_c1
+  nrm_len = wp.length(diff)
+  if nrm_len > 1e-10:
+    normal = diff / nrm_len
+  else:
+    e1 = a2 - a1
+    e2 = a3 - a1
+    normal = wp.normalize(wp.cross(e1, e2))
+
+  pos = best_c1 + normal * (a_radius + 0.5 * dist)
+  return dist, pos, normal
+
+
 @wp.func
 def box_triangle(
   box_pos: wp.vec3,
diff --git a/mujoco_warp/_src/collision_primitive_core_test.py b/mujoco_warp/_src/collision_primitive_core_test.py
index 206670566..55f220b79 100644
--- a/mujoco_warp/_src/collision_primitive_core_test.py
+++ b/mujoco_warp/_src/collision_primitive_core_test.py
@@ -531,5 +531,107 @@ def test_with_triangle_radius(self):
     self.assertLess(dist[0], collision_primitive_core.MJ_MAXVAL)
 
 
+@wp.kernel
+def triangle_triangle_kernel(
+  # In:
+  a1: wp.vec3,
+  a2: wp.vec3,
+  a3: wp.vec3,
+  a_radius: float,
+  b1: wp.vec3,
+  b2: wp.vec3,
+  b3: wp.vec3,
+  b_radius: float,
+  # Out:
+  dist_out: wp.array(dtype=float),
+  pos_out: wp.array(dtype=wp.vec3),
+  normal_out: wp.array(dtype=wp.vec3),
+):
+  dist, pos, normal = collision_primitive_core.triangle_triangle(a1, a2, a3, a_radius, b1, b2, b3, b_radius)
+  dist_out[0] = dist
+  pos_out[0] = pos
+  normal_out[0] = normal
+
+
+class TriangleTriangleTest(parameterized.TestCase):
+  def _run(self, a1, a2, a3, a_r, b1, b2, b3, b_r):
+    dist = wp.zeros(1, dtype=float)
+    pos = wp.zeros(1, dtype=wp.vec3)
+    normal = wp.zeros(1, dtype=wp.vec3)
+    wp.launch(
+      triangle_triangle_kernel,
+      dim=1,
+      inputs=[
+        wp.vec3(a1),
+        wp.vec3(a2),
+        wp.vec3(a3),
+        a_r,
+        wp.vec3(b1),
+        wp.vec3(b2),
+        wp.vec3(b3),
+        b_r,
+      ],
+      outputs=[dist, pos, normal],
+    )
+    return dist.numpy()[0], pos.numpy()[0], normal.numpy()[0]
+
+  def test_separated_parallel(self):
+    a1 = np.array([0.0, 0.0, 0.0])
+    a2 = np.array([1.0, 0.0, 0.0])
+    a3 = np.array([0.5, 1.0, 0.0])
+    b1 = np.array([0.0, 0.0, 2.0])
+    b2 = np.array([1.0, 0.0, 2.0])
+    b3 = np.array([0.5, 1.0, 2.0])
+    dist, _, _ = self._run(a1, a2, a3, 0.0, b1, b2, b3, 0.0)
+    self.assertGreater(dist, 1.0)
+
+  def test_vertex_near_face(self):
+    a1 = np.array([0.0, 0.0, 0.0])
+    a2 = np.array([1.0, 0.0, 0.0])
+    a3 = np.array([0.5, 1.0, 0.0])
+    b1 = np.array([0.4, 0.3, 0.1])
+    b2 = np.array([0.4, 0.3, 1.0])
+    b3 = np.array([1.0, 0.3, 1.0])
+    r = 0.02
+    dist, _, normal = self._run(a1, a2, a3, r, b1, b2, b3, r)
+    expected_dist = 0.1 - 2 * r
+    np.testing.assert_allclose(dist, expected_dist, atol=1e-4)
+    self.assertAlmostEqual(abs(normal[2]), 1.0, places=3)
+
+  def test_penetrating(self):
+    a1 = np.array([0.0, 0.0, 0.0])
+    a2 = np.array([1.0, 0.0, 0.0])
+    a3 = np.array([0.5, 1.0, 0.0])
+    b1 = np.array([0.4, 0.3, -0.03])
+    b2 = np.array([0.4, 0.3, 1.0])
+    b3 = np.array([1.0, 0.3, 1.0])
+    r = 0.02
+    dist, _, _ = self._run(a1, a2, a3, r, b1, b2, b3, r)
+    self.assertLess(dist, 0.0)
+
+  def test_with_radius_closes_gap(self):
+    a1 = np.array([0.0, 0.0, 0.0])
+    a2 = np.array([1.0, 0.0, 0.0])
+    a3 = np.array([0.5, 1.0, 0.0])
+    b1 = np.array([0.4, 0.3, 0.08])
+    b2 = np.array([0.4, 0.3, 1.0])
+    b3 = np.array([1.0, 0.3, 1.0])
+    dist_no_r, _, _ = self._run(a1, a2, a3, 0.0, b1, b2, b3, 0.0)
+    dist_with_r, _, _ = self._run(a1, a2, a3, 0.05, b1, b2, b3, 0.05)
+    self.assertGreater(dist_no_r, 0.0)
+    self.assertLess(dist_with_r, 0.0)
+
+  def test_edge_proximity(self):
+    a1 = np.array([0.0, 0.0, 0.0])
+    a2 = np.array([1.0, 0.0, 0.0])
+    a3 = np.array([0.5, 1.0, 0.0])
+    b1 = np.array([0.5, -0.1, 0.15])
+    b2 = np.array([0.5, -0.1, 1.0])
+    b3 = np.array([1.5, -0.1, 0.5])
+    r = 0.02
+    dist, _, _ = self._run(a1, a2, a3, r, b1, b2, b3, r)
+    self.assertLess(dist, collision_primitive_core.MJ_MAXVAL)
+
+
 if __name__ == "__main__":
   absltest.main()
diff --git a/mujoco_warp/_src/io.py b/mujoco_warp/_src/io.py
index a923917b3..3f5f02e92 100644
--- a/mujoco_warp/_src/io.py
+++ b/mujoco_warp/_src/io.py
@@ -223,6 +223,23 @@ def _check_friction(name: str, id_: int, condim: int, friction, checks):
   m.is_sparse = is_sparse(mjm)
   m.has_fluid = mjm.opt.wind.any() or mjm.opt.density > 0 or mjm.opt.viscosity > 0
 
+  if mjm.nflex > 0:
+    selfcollide = np.array(mjm.flex_selfcollide)
+    unsupported = (selfcollide != 0) & (selfcollide != 1)
+    if unsupported.any():
+      raise NotImplementedError(
+        'flex_selfcollide modes BVH/SAP/AUTO not supported in MuJoCo-Warp. Use selfcollide="narrow" or selfcollide="none".'
+      )
+    # Precompute max self-collision pairs for dim=2 elements to avoid
+    # device→host copies during CUDA graph capture.
+    n_dim2_elems = 0
+    for f in range(mjm.nflex):
+      if selfcollide[f] != 0 and mjm.flex_dim[f] == 2:
+        n_dim2_elems += mjm.flex_elemnum[f]
+    m.flex_self_max_pairs = n_dim2_elems * (n_dim2_elems - 1) // 2
+  else:
+    m.flex_self_max_pairs = 0
+
   m.max_ten_J_rownnz = int(mjm.ten_J_rownnz.max()) if mjm.ntendon else 0
 
   # body ids grouped by tree level (depth-based traversal)
@@ -1087,6 +1104,11 @@ def put_data(
     else:
       njmax_nnz = njmax * mjm.nv
 
+  # ensure njmax_nnz can hold the actual data from mjd
+  if is_sparse(mjm) and mjd.nefc > 0 and mujoco.mj_isSparse(mjm):
+    actual_nnz = int(mjd.efc_J_rownnz[: mjd.nefc].sum())
+    njmax_nnz = max(njmax_nnz, actual_nnz)
+
   # ensure static geom positions are computed
   # TODO: remove once MjData creation semantics are fixed
   mujoco.mj_kinematics(mjm, mjd)
diff --git a/mujoco_warp/_src/smooth.py b/mujoco_warp/_src/smooth.py
index 7354e16c6..866a618ea 100644
--- a/mujoco_warp/_src/smooth.py
+++ b/mujoco_warp/_src/smooth.py
@@ -354,6 +354,65 @@ def _flex_edges(
       flexedge_J_out[worldid, rowadr + nnz_offset + k] = wp.dot(jacp2, edge)
 
 
+@wp.kernel
+def _compute_flexelem_aabb(
+  # Model:
+  nflex: int,
+  flex_dim: wp.array(dtype=int),
+  flex_vertadr: wp.array(dtype=int),
+  flex_elemadr: wp.array(dtype=int),
+  flex_elemnum: wp.array(dtype=int),
+  flex_elem: wp.array(dtype=int),
+  flex_radius: wp.array(dtype=float),
+  # Data in:
+  flexvert_xpos_in: wp.array2d(dtype=wp.vec3),
+  # Data out:
+  flexelem_aabb_out: wp.array3d(dtype=float),
+):
+  worldid, elemid = wp.tid()
+
+  flexid = int(-1)
+  for i in range(nflex):
+    elem_adr = flex_elemadr[i]
+    elem_num = flex_elemnum[i]
+    if elemid >= elem_adr and elemid < elem_adr + elem_num:
+      flexid = i
+      break
+
+  if flexid < 0:
+    return
+
+  dim = flex_dim[flexid]
+  vertadr = flex_vertadr[flexid]
+  radius = flex_radius[flexid]
+
+  # Compute element data offset by accumulating elem counts for preceding
+  # flexes, since elements of different dims are packed contiguously in
+  # flex_elem with strides of (dim+1).
+  elem_data_base = int(0)
+  for i in range(flexid):
+    elem_data_base += flex_elemnum[i] * (flex_dim[i] + 1)
+  local_elem = elemid - flex_elemadr[flexid]
+  elem_data_idx = elem_data_base + local_elem * (dim + 1)
+
+  v0 = flex_elem[elem_data_idx]
+
+  v0 = flex_elem[elemid * (dim + 1)]
+  p0 = flexvert_xpos_in[worldid, vertadr + v0]
+  xmin = wp.vec3(p0[0], p0[1], p0[2])
+  xmax = wp.vec3(p0[0], p0[1], p0[2])
+
+  for i in range(1, dim + 1):
+    vi = flex_elem[elem_data_idx + i]
+    pi = flexvert_xpos_in[worldid, vertadr + vi]
+    xmin = wp.vec3(wp.min(xmin[0], pi[0]), wp.min(xmin[1], pi[1]), wp.min(xmin[2], pi[2]))
+    xmax = wp.vec3(wp.max(xmax[0], pi[0]), wp.max(xmax[1], pi[1]), wp.max(xmax[2], pi[2]))
+
+  for axis in range(3):
+    flexelem_aabb_out[worldid, elemid, axis] = 0.5 * (xmax[axis] + xmin[axis])
+    flexelem_aabb_out[worldid, elemid, axis + 3] = 0.5 * (xmax[axis] - xmin[axis]) + radius
+
+
 @event_scope
 def kinematics(m: Model, d: Data):
   """Computes forward kinematics for all bodies, sites, geoms, and flexible elements.
@@ -459,6 +518,25 @@ def flex(m: Model, d: Data):
     ],
   )
 
+  if m.nflexelem > 0:
+    wp.launch(
+      _compute_flexelem_aabb,
+      dim=(d.nworld, m.nflexelem),
+      inputs=[
+        m.nflex,
+        m.flex_dim,
+        m.flex_vertadr,
+        m.flex_elemadr,
+        m.flex_elemnum,
+        m.flex_elem,
+        m.flex_radius,
+        d.flexvert_xpos,
+      ],
+      outputs=[
+        d.flexelem_aabb,
+      ],
+    )
+
 
 @wp.kernel
 def _subtree_com_init(
diff --git a/mujoco_warp/_src/smooth_test.py b/mujoco_warp/_src/smooth_test.py
index 92dfa1ab6..25cad23d1 100644
--- a/mujoco_warp/_src/smooth_test.py
+++ b/mujoco_warp/_src/smooth_test.py
@@ -563,7 +563,7 @@ def test_flex_1d_pinned(self):
           <geom type="sphere" size="0.02" mass="0.01"/>
           <flexcomp name="line" type="grid" count="5 1 1" spacing="0.1 0.1 0.1"
                     radius="0.01" dim="1" mass="1">
-            <contact contype="0" conaffinity="0"/>
+            <contact contype="0" conaffinity="0" selfcollide="none"/>
             <edge equality="true" damping="0.01"/>
             <pin id="0"/>
           </flexcomp>
diff --git a/mujoco_warp/_src/types.py b/mujoco_warp/_src/types.py
index 174ec3f65..9495dc849 100644
--- a/mujoco_warp/_src/types.py
+++ b/mujoco_warp/_src/types.py
@@ -96,6 +96,24 @@ class BroadphaseFilter(enum.IntFlag):
   OBB = 8
 
 
+class FlexSelfCollideType(enum.IntEnum):
+  """Type of flex self-collision algorithm.
+
+  Attributes:
+    NONE: No self-collisions
+    NARROW: Skip broadphase, go directly to narrowphase (brute force)
+    BVH: Use BVH in midphase
+    SAP: Use sweep-and-prune in midphase
+    AUTO: Choose automatically based on flex dimension
+  """
+
+  NONE = 0
+  NARROW = 1
+  BVH = 2
+  SAP = 3
+  AUTO = 4
+
+
 class CamLightType(enum.IntEnum):
   """Type of camera light.
 
@@ -973,8 +991,11 @@ class Model:
     flex_contype: flex contact type                          (nflex,)
     flex_conaffinity: flex contact affinity                  (nflex,)
     flex_condim: contact dimensionality (1, 3, 4, 6)         (nflex,)
+    flex_solref: solver reference                            (nflex, 2)
+    flex_solimp: solver impedance                            (nflex, 5)
     flex_friction: friction for (slide, spin, roll)          (nflex, 3)
     flex_margin: detect contact if dist<margin               (nflex,)
+    flex_selfcollide: self-collision mode                    (nflex,)
     flex_dim: 1: lines, 2: triangles, 3: tetrahedra          (nflex,)
     flex_vertadr: first vertex address                       (nflex,)
     flex_vertnum: number of vertices                         (nflex,)
@@ -1150,6 +1171,7 @@ class Model:
     geom_pair_type_count: count of max number of each
                           potential collision
     geom_plugin_index: geom index in plugin array            (ngeom,)
+    flex_self_max_pairs: max number of flex self contact pairs
     eq_connect_adr: eq_* addresses of type `CONNECT`
     eq_wld_adr: eq_* addresses of type `WELD`
     eq_jnt_adr: eq_* addresses of type `JOINT`
@@ -1360,8 +1382,11 @@ class Model:
   flex_contype: array("nflex", int)
   flex_conaffinity: array("nflex", int)
   flex_condim: array("nflex", int)
+  flex_solref: array("nflex", wp.vec2)
+  flex_solimp: array("nflex", vec5)
   flex_friction: array("nflex", wp.vec3)
   flex_margin: array("nflex", float)
+  flex_selfcollide: array("nflex", int)
   flex_dim: array("nflex", int)
   flex_vertadr: array("nflex", int)
   flex_vertnum: array("nflex", int)
@@ -1531,6 +1556,7 @@ class Model:
   nxn_pairid_filtered: wp.array(dtype=wp.vec2i)
   geom_pair_type_count: tuple[int, ...]
   geom_plugin_index: array("ngeom", int)
+  flex_self_max_pairs: int
   eq_connect_adr: wp.array(dtype=int)
   eq_wld_adr: wp.array(dtype=int)
   eq_jnt_adr: wp.array(dtype=int)
@@ -1721,6 +1747,7 @@ class Data:
     cdof: com-based motion axis of each dof (rot:lin)           (nworld, nv, 6)
     cinert: com-based body inertia and mass                     (nworld, nbody, 10)
     flexvert_xpos: cartesian flex vertex positions              (nworld, nflexvert, 3)
+    flexelem_aabb: element bounding boxes (center, half-size)   (nworld, nflexelem, 6)
     flexedge_J: edge length Jacobian                            (nworld, nJfe)
     flexedge_length: flex edge lengths                          (nworld, nflexedge, 1)
     ten_wrapadr: start address of tendon's path                 (nworld, ntendon)
@@ -1819,6 +1846,7 @@ class Data:
   cdof: array("nworld", "nv", wp.spatial_vector)
   cinert: array("nworld", "nbody", vec10)
   flexvert_xpos: array("nworld", "nflexvert", wp.vec3)
+  flexelem_aabb: array("nworld", "nflexelem", 6, float)
   flexedge_J: array("nworld", "nJfe", float)
   flexedge_length: array("nworld", "nflexedge", float)
   ten_wrapadr: array("nworld", "ntendon", int)