Scale flex warp constraint by thickness^3.

Warp is not a true stretch mode but a parasitic shear mode so it should be scaled like bending (~thickness^3) rather than stretching (~thickness). The thickness parameter can therefore be used to set the stretch to warp stiffness ratio of the strain constraints.

PiperOrigin-RevId: 916900311
Change-Id: I52e54210fbb2a6d309abc1f084934830adce3423

model/flex/hollow_vs_solid.xml

@@ -50,7 +50,7 @@
         origin="0 0 0" count="8 2 12" cellcount="6 1 6"
         type="grid" name="soft_mesh_2" dim="3" spacing=".025 .05 .025"
         mass="0.43" radius="0.005" dof="trilinear" rgba="0.9 0.7 0.7 1">
-        <elasticity young="1e4" poisson="0.3" damping="0" elastic2d="bend" thickness="0.03"/>
+        <elasticity young="1e4" poisson="0.3" damping="0" elastic2d="bend" thickness="0.02"/>
         <edge equality="strain" solref="0.001 1" solimp="0.99 0.999 0.001"/>
         <contact selfcollide="none"
               internal="false"

src/user/user_mesh.cc

@@ -4071,8 +4071,9 @@ static void ComputeWarpMode(double* warp, const double* pos,
 }
 
 
-// compute the warp stiffness for a 2D face element, matching the transverse
+// compute the warp bending stiffness for a 2D face element
-// shear energy of the old 3D-on-2D formulation: λ_warp = μ t (d1/d0 + d0/d1)/6
+// uses plate bending theory: the warp mode is a pure twist (κ_xy),
+// with bending stiffness proportional to t³ (no shear locking)
 //   pos:          node positions (3*npe doubles)
 //   npe:          nodes per element
 //   normal_axis:  axis perpendicular to the face
@@ -4082,13 +4083,17 @@ static double ComputeWarpStiffness(const double* pos, int npe, int normal_axis,
                                    double E, double nu, double thickness) {
   int axis0 = (normal_axis + 1) % 3;
   int axis1 = (normal_axis + 2) % 3;
-  double d0 = pos[3*(npe-1) + axis0] - pos[axis0];
+  double d0 = std::abs(pos[3*(npe-1) + axis0] - pos[axis0]);
-  double d1 = pos[3*(npe-1) + axis1] - pos[axis1];
+  double d1 = std::abs(pos[3*(npe-1) + axis1] - pos[axis1]);
-  double mu = E / (2.0 * (1.0 + nu));
 
-  // warp stiffness from transverse shear Rayleigh quotient:
+  if (d0 < 1e-30 || d1 < 1e-30) return 0;
-  //   w^T K_phys w / |w|^2 = μ t (|d1/d0| + |d0/d1|) / 6
+
-  return mu * thickness * (std::abs(d1/d0) + std::abs(d0/d1)) / 6.0;
+  // plate bending rigidity: D = E*t³ / (12*(1-ν²))
+  double D = E * thickness * thickness * thickness / (12.0 * (1.0 - nu * nu));
+
+  // warp stiffness from twist curvature Rayleigh quotient:
+  //   w^T K_bend w / |w|^2 = D*(1-ν)*4 / (d0*d1)
+  return D * (1.0 - nu) * 4.0 / (d0 * d1);
 }
 
 
@@ -5069,10 +5074,10 @@ void mjCFlex::Compile(const mjVFS* vfs) {
 
         if (shell_mode) {
           // pure membrane K: eigendecompose gives 5 membrane modes (Q1),
-          // then we add 1 explicit warp mode
+          // then we add 1 explicit warp mode with bending stiffness (∝ t³)
           int neig = EigendecomposeStiffness(K_elem.data(), out, ndof_elem);
 
-          // add explicit warp mode with stiffness matching old transverse shear
+          // add explicit warp mode with plate bending stiffness
           double warp_stiffness = ComputeWarpStiffness(
               elem_pos.data(), npe, normal_axis, K_young, K_poisson, thickness);
           if (warp_stiffness > 0) {