Bending stiffness for curved shells in flex.

quagla · copybara-github · commit b66175eba62c · 2025-08-18T04:13:18.000-07:00
PiperOrigin-RevId: 796361841
Change-Id: I105ea235fe2a8e1bcbde67276d77fa92eed214a1
diff --git a/doc/changelog.rst b/doc/changelog.rst
@@ -5,6 +5,11 @@ Changelog
 Upcoming version (not yet released)
 -----------------------------------
 
+General
+^^^^^^^
+- Added support for shells with a curved reference configuration. See this `example
+   <https://github.com/google-deepmind/mujoco/blob/main/model/flex/basket.xml>`__
+
 MJX
 ^^^
 - Promote ``ten_length`` to the public MJX API. Add Warp support for ``mjx.tendon``.
diff --git a/doc/includes/references.h b/doc/includes/references.h
@@ -1268,7 +1268,7 @@ struct mjModel_ {
   mjtNum*   flexedge_invweight0;  // edge inv. weight in qpos0                (nflexedge x 1)
   mjtNum*   flex_radius;          // radius around primitive element          (nflex x 1)
   mjtNum*   flex_stiffness;       // finite element stiffness matrix          (nflexelem x 21)
-  mjtNum*   flex_bending;         // bending stiffness                        (nflexedge x 16)
+  mjtNum*   flex_bending;         // bending stiffness                        (nflexedge x 17)
   mjtNum*   flex_damping;         // Rayleigh's damping coefficient           (nflex x 1)
   mjtNum*   flex_edgestiffness;   // edge stiffness                           (nflex x 1)
   mjtNum*   flex_edgedamping;     // edge damping                             (nflex x 1)
diff --git a/include/mujoco/mjmodel.h b/include/mujoco/mjmodel.h
@@ -958,7 +958,7 @@ struct mjModel_ {
   mjtNum*   flexedge_invweight0;  // edge inv. weight in qpos0                (nflexedge x 1)
   mjtNum*   flex_radius;          // radius around primitive element          (nflex x 1)
   mjtNum*   flex_stiffness;       // finite element stiffness matrix          (nflexelem x 21)
-  mjtNum*   flex_bending;         // bending stiffness                        (nflexedge x 16)
+  mjtNum*   flex_bending;         // bending stiffness                        (nflexedge x 17)
   mjtNum*   flex_damping;         // Rayleigh's damping coefficient           (nflex x 1)
   mjtNum*   flex_edgestiffness;   // edge stiffness                           (nflex x 1)
   mjtNum*   flex_edgedamping;     // edge damping                             (nflex x 1)
diff --git a/include/mujoco/mjxmacro.h b/include/mujoco/mjxmacro.h
@@ -372,7 +372,7 @@
     X   ( mjtNum,  flexedge_invweight0,   nflexedge,     1                    ) \
     X   ( mjtNum,  flex_radius,           nflex,         1                    ) \
     X   ( mjtNum,  flex_stiffness,        nflexelem,     21                   ) \
-    X   ( mjtNum,  flex_bending,          nflexedge,     16                   ) \
+    X   ( mjtNum,  flex_bending,          nflexedge,     17                   ) \
     X   ( mjtNum,  flex_damping,          nflex,         1                    ) \
     X   ( mjtNum,  flex_edgestiffness,    nflex,         1                    ) \
     X   ( mjtNum,  flex_edgedamping,      nflex,         1                    ) \
diff --git a/mjx/mujoco/mjx/third_party/mujoco_warp/_src/io.py b/mjx/mujoco/mjx/third_party/mujoco_warp/_src/io.py
@@ -609,7 +609,7 @@ def create_nmodel_batched_array(mjm_array, dtype, expand_dim=True):
     flex_elemedge=wp.array(mjm.flex_elemedge, dtype=int),
     flexedge_length0=wp.array(mjm.flexedge_length0, dtype=float),
     flex_stiffness=wp.array(mjm.flex_stiffness.flatten(), dtype=float),
-    flex_bending=wp.array(mjm.flex_bending, dtype=wp.mat44f),
+    flex_bending=wp.array(mjm.flex_bending.flatten(), dtype=float),
     flex_damping=wp.array(mjm.flex_damping, dtype=float),
     mesh_vertadr=wp.array(mjm.mesh_vertadr, dtype=int),
     mesh_vertnum=wp.array(mjm.mesh_vertnum, dtype=int),
diff --git a/mjx/mujoco/mjx/third_party/mujoco_warp/_src/passive.py b/mjx/mujoco/mjx/third_party/mujoco_warp/_src/passive.py
@@ -454,7 +454,7 @@ def _flex_bending(
   flex_vertbodyid: wp.array(dtype=int),
   flex_edge: wp.array(dtype=wp.vec2i),
   flex_edgeflap: wp.array(dtype=wp.vec2i),
-  flex_bending: wp.array(dtype=wp.mat44f),
+  flex_bending: wp.array(dtype=float),
   # Data in:
   flexvert_xpos_in: wp.array2d(dtype=wp.vec3),
   # Data out:
@@ -481,7 +481,7 @@ def _flex_bending(
   for i in range(nvert):
     for j in range(nvert):
       for x in range(3):
-        force[i, x] -= flex_bending[edgeid][i, j] * flexvert_xpos_in[worldid, v[j]][x]
+        force[i, x] -= flex_bending[17*edgeid + 4*i + j] * flexvert_xpos_in[worldid, v[j]][x]
 
   for i in range(nvert):
     bodyid = flex_vertbodyid[flex_vertadr[f] + v[i]]
diff --git a/mjx/mujoco/mjx/third_party/mujoco_warp/_src/types.py b/mjx/mujoco/mjx/third_party/mujoco_warp/_src/types.py
@@ -1199,7 +1199,7 @@ class Model:
   flex_elemedge: wp.array(dtype=int)
   flexedge_length0: wp.array(dtype=float)
   flex_stiffness: wp.array(dtype=float)
-  flex_bending: wp.array(dtype=wp.mat44f)
+  flex_bending: wp.array(dtype=float)
   flex_damping: wp.array(dtype=float)
   mesh_vertadr: wp.array(dtype=int)
   mesh_vertnum: wp.array(dtype=int)
diff --git a/mjx/mujoco/mjx/warp/forward.py b/mjx/mujoco/mjx/warp/forward.py
@@ -124,7 +124,7 @@ def _forward_shim(
     eq_solref: wp.array2d(dtype=wp.vec2),
     eq_ten_adr: wp.array(dtype=int),
     eq_wld_adr: wp.array(dtype=int),
-    flex_bending: wp.array(dtype=wp.mat44f),
+    flex_bending: wp.array(dtype=float),
     flex_damping: wp.array(dtype=float),
     flex_dim: wp.array(dtype=int),
     flex_edge: wp.array(dtype=wp.vec2i),
@@ -2148,7 +2148,7 @@ def _step_shim(
     eq_solref: wp.array2d(dtype=wp.vec2),
     eq_ten_adr: wp.array(dtype=int),
     eq_wld_adr: wp.array(dtype=int),
-    flex_bending: wp.array(dtype=wp.mat44f),
+    flex_bending: wp.array(dtype=float),
     flex_damping: wp.array(dtype=float),
     flex_dim: wp.array(dtype=int),
     flex_edge: wp.array(dtype=wp.vec2i),
diff --git a/mjx/mujoco/mjx/warp/types.py b/mjx/mujoco/mjx/warp/types.py
@@ -737,7 +737,7 @@ def _from_elt(cont, axis_size, d, axis_dest):
         'eq_type': 1,
         'eq_wld_adr': 1,
         'exclude_signature': 1,
-        'flex_bending': 3,
+        'flex_bending': 1,
         'flex_damping': 1,
         'flex_dim': 1,
         'flex_edge': 2,
diff --git a/model/flex/basket.xml b/model/flex/basket.xml
@@ -0,0 +1,48 @@
+<!-- Copyright 2025 DeepMind Technologies Limited
+
+     Licensed under the Apache License, Version 2.0 (the "License");
+     you may not use this file except in compliance with the License.
+     You may obtain a copy of the License at
+
+         http://www.apache.org/licenses/LICENSE-2.0
+
+     Unless required by applicable law or agreed to in writing, software
+     distributed under the License is distributed on an "AS IS" BASIS,
+     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+     See the License for the specific language governing permissions and
+     limitations under the License.
+-->
+
+<mujoco model="Basket">
+  <include file="scene.xml"/>
+
+  <compiler autolimits="true"/>
+
+  <option solver="CG" tolerance="1e-6"/>
+
+  <size memory="10M"/>
+
+  <visual>
+    <map stiffness="100"/>
+  </visual>
+
+  <worldbody>
+    <flexcomp type="box" count="8 8 8" spacing=".1 .1 .1" pos="0 0 1"
+              radius=".01" rgba=".68 .53 .38 1" name="box" dim="2" mass="0.25">
+      <contact condim="3" solref="0.01 1" solimp=".95 .99 .0001" selfcollide="none"/>
+      <edge equality="true" damping="10"/>
+      <elasticity young="6e6" poisson="0.2" thickness="8e-3" elastic2d="bend"/>
+    </flexcomp>
+
+    <replicate count="3" offset=".2 0 0">
+      <replicate count="3" offset="0 .2 0">
+        <replicate count="5" offset="0 0 .2">
+          <body name="ball" pos="-.2 -.2 1.8">
+            <freejoint/>
+            <geom type="sphere" size="0.1" rgba="1 0 .2 1" mass=".03"/>
+          </body>
+        </replicate>
+      </replicate>
+    </replicate>
+  </worldbody>
+</mujoco>
diff --git a/python/mujoco/introspect/structs.py b/python/mujoco/introspect/structs.py
@@ -2811,7 +2811,7 @@
                      inner_type=ValueType(name='mjtNum'),
                  ),
                  doc='bending stiffness',
-                 array_extent=('nflexedge', 16),
+                 array_extent=('nflexedge', 17),
              ),
              StructFieldDecl(
                  name='flex_damping',
diff --git a/src/engine/engine_passive.c b/src/engine/engine_passive.c
@@ -116,7 +116,7 @@ static void mj_springdamper(const mjModel* m, mjData* d) {
   // flex elasticity
   for (int f=0; f < m->nflex; f++) {
     mjtNum* k = m->flex_stiffness + 21*m->flex_elemadr[f];
-    mjtNum* b = m->flex_bending + 16*m->flex_edgeadr[f];
+    mjtNum* b = m->flex_bending + 17*m->flex_edgeadr[f];
     int dim = m->flex_dim[f];
 
     if (dim == 1 || m->flex_rigid[f]) {
@@ -136,12 +136,32 @@ static void mj_springdamper(const mjModel* m, mjData* d) {
           // skip boundary edges
           continue;
         }
+
+        // flap edges
+        mjtNum ed[3][3];
+        mju_sub3(ed[0], xpos + 3*v[1], xpos + 3*v[0]);
+        mju_sub3(ed[1], xpos + 3*v[2], xpos + 3*v[0]);
+        mju_sub3(ed[2], xpos + 3*v[3], xpos + 3*v[0]);
+
+        // forces at the vertices due to curved reference
+        mjtNum frc[4][3];
+        mju_cross(frc[1], ed[1], ed[2]);
+        mju_cross(frc[2], ed[2], ed[0]);
+        mju_cross(frc[3], ed[0], ed[1]);
+        frc[0][0] = -(frc[1][0] + frc[2][0] + frc[3][0]);
+        frc[0][1] = -(frc[1][1] + frc[2][1] + frc[3][1]);
+        frc[0][2] = -(frc[1][2] + frc[2][2] + frc[3][2]);
+
+        // force
         mjtNum force[12] = {0};
-        for (int i = 0; i < 4; i++) {
-          for (int j = 0; j < 4; j++) {
-            for (int x = 0; x < 3; x++) {
-              force[3*i+x] += b[16*e+4*i+j] * xpos[3*v[j]+x];
+        for (int x = 0; x < 3; x++) {
+          for (int i = 0; i < 4; i++) {
+            for (int j = 0; j < 4; j++) {
+              // thin plate bending force
+              force[3*i+x] += b[17*e+4*i+j] * xpos[3*v[j]+x];
             }
+            // curved reference contribution
+            force[3*i+x] += b[17*e+16] * frc[i][x];
           }
         }
 
diff --git a/src/user/user_flexcomp.cc b/src/user/user_flexcomp.cc
diff --git a/src/user/user_flexcomp.h b/src/user/user_flexcomp.h
diff --git a/src/user/user_mesh.cc b/src/user/user_mesh.cc
diff --git a/src/user/user_model.cc b/src/user/user_model.cc
diff --git a/test/engine/engine_passive_test.cc b/test/engine/engine_passive_test.cc
diff --git a/test/plugin/elasticity/elasticity_test.cc b/test/plugin/elasticity/elasticity_test.cc
