recommended stable VBD cloth parameters and grasp logic for towel diagonal folding

[N3rv3rM1nd]

Hi Newton team and community members,

Sorry for disturbing you.

First of all, thank you very much for open-sourcing Newton and for providing the cloth examples. I am currently using Newton for my project on deformable object manipulation, and I would sincerely appreciate some guidance.

My target task is diagonal folding of a flat towel. In the long run I hope to train a policy for this task, but at the current stage I am still trying to build a stable simulation setup for the cloth itself.

I have already tried to tune the cloth setup by referring to the official Newton cloth examples and their parameter styles. I also looked at the example that seems most relevant to my task, cloth_franka. However, I feel that this example is still somewhat different from my scenario: it uses a garment-like shape rather than a flat towel, and it is set up in centimeter scale, so I have not yet found a good way to transfer that setup to my current towel-folding task.

At the moment, my main difficulty is that I have not yet found a stable combination of:

1. towel-like cloth parameters in VBD, and
2. a suitable grasp logic for lifting and folding the cloth.

I have already tried several directions on my side:

• tuning cloth material parameters such as density, tri_ke, tri_ka, tri_kd, edge_ke, edge_kd, particle radius, resolution, and self-contact settings;
• simplifying the problem into cloth-only branch experiments, without robot dynamics;
• trying different virtual grasp strategies, including direct particle/vertex driving, small patch-based grasping, and contact-only kinematic pinch style grasping.

However, I am still facing two opposite difficulties:

• if I directly drive cloth particles or a small patch, the cloth often becomes numerically unstable during lifting, with severe local distortion and unrealistic deformation;
• if I rely only on contact-based pinch style grasping, the cloth becomes much more stable, but it often cannot be lifted effectively.

Because of this, I would like to ask:

• Is there a recommended parameter regime in Newton VBD for towel-like cloth?
• Is there a grasp logic that you would recommend for cloth manipulation / folding tasks in Newton?
• In your experience, is direct particle driving generally too aggressive for this kind of task?
• Is a hybrid approach (for example, contact-gated virtual grasp or weak attachment logic) a more appropriate direction?

I am not asking for a complete solution, but I would be very grateful for any practical guidance on what Newton currently considers a stable setup for cloth folding or towel manipulation.

If it would be helpful, I would also be happy to prepare and share a minimal reproducible cloth-only script.

Thank you very much for your time and help.

Codes are below:

import math
import numpy as np
import warp as wp
import newton
import newton.examples
from newton.solvers import SolverVBD

wp.init()


def build_unique_edges_from_triangles(indices: list[int]) -> np.ndarray:
    edges = set()
    for t in range(0, len(indices), 3):
        a, b, c = int(indices[t]), int(indices[t+1]), int(indices[t+2])
        for i, j in ((a, b), (b, c), (c, a)):
            if i > j:
                i, j = j, i
            edges.add((i, j))
    return np.array(sorted(edges), dtype=np.int32)


def stretch_stats(particles: np.ndarray, edges: np.ndarray, rest_lengths: np.ndarray):
    if len(edges) == 0:
        return float('nan'), float('nan'), float('nan'), np.array([], dtype=np.float64)
    cur = np.linalg.norm(particles[edges[:, 0]] - particles[edges[:, 1]], axis=1)
    ratio = cur / np.maximum(rest_lengths, 1e-12)
    return float(np.max(ratio)), float(np.percentile(ratio, 95.0)), float(np.mean(ratio)), ratio


def build_triangles(indices: list[int]) -> np.ndarray:
    return np.asarray(indices, dtype=np.int32).reshape(-1, 3)


def triangle_rest_geometry(rest_particles: np.ndarray, tris: np.ndarray):
    a = rest_particles[tris[:, 0]]
    b = rest_particles[tris[:, 1]]
    c = rest_particles[tris[:, 2]]
    rest_cross = np.cross(b - a, c - a)
    rest_area = 0.5 * np.linalg.norm(rest_cross, axis=1)
    rest_normal_unit = rest_cross / np.maximum(np.linalg.norm(rest_cross, axis=1, keepdims=True), 1e-12)
    return rest_area, rest_normal_unit


def triangle_signed_area_stats(particles: np.ndarray, tris: np.ndarray, rest_area: np.ndarray, rest_normal_unit: np.ndarray):
    if len(tris) == 0:
        return (float('nan'), float('nan'), float('nan'), 0, np.array([], dtype=np.float64))
    a = particles[tris[:, 0]]
    b = particles[tris[:, 1]]
    c = particles[tris[:, 2]]
    cur_cross = np.cross(b - a, c - a)
    signed_double_area = np.sum(cur_cross * rest_normal_unit, axis=1)
    signed_ratio = signed_double_area / np.maximum(2.0 * rest_area, 1e-12)
    flip_count = int(np.sum(signed_ratio < 0.0))
    return (
        float(np.min(signed_ratio)),
        float(np.percentile(signed_ratio, 5.0)),
        float(np.mean(signed_ratio)),
        flip_count,
        signed_ratio,
    )


def build_corner_topology_patch(corner_idx: int, res: int, side: int = 2) -> np.ndarray:
    """Create a fixed topological corner patch (e.g. 2x2 or 3x3) that extends inward from the detected corner."""
    ii, jj = divmod(int(corner_idx), res)
    step_i = 1 if ii < (res // 2) else -1
    step_j = 1 if jj < (res // 2) else -1

    patch = []
    for di in range(side):
        for dj in range(side):
            i = ii + di * step_i
            j = jj + dj * step_j
            i = max(0, min(res - 1, i))
            j = max(0, min(res - 1, j))
            patch.append(i * res + j)
    return np.array(sorted(set(patch)), dtype=np.int32)


def select_target_patch_spatial(particles: np.ndarray, query_xy: np.ndarray, radius_xy=0.012, n_pts=6):
    xy = particles[:, :2]
    dxy = np.linalg.norm(xy - query_xy[None, :], axis=1)
    cand = np.where(dxy < radius_xy)[0]
    if len(cand) < n_pts:
        cand = np.argsort(dxy)[:n_pts]
    else:
        cand = cand[np.argsort(dxy[cand])[:n_pts]]
    return np.array(cand, dtype=np.int32)


class Example:
    def __init__(self, viewer, args):
        self.viewer = viewer
        self.args = args
        self.sim_time = 0.0
        self.fps = 60
        self.frame_dt = 1.0 / self.fps
        self.sim_substeps = 20
        self.sim_dt = self.frame_dt / self.sim_substeps

        # Cloth / scene
        self.cloth_size = 0.15
        self.cloth_res = 31
        self.table_center = wp.vec3(0.6, 0.0, 0.05)
        self.table_half = (0.4, 0.5, 0.05)
        self.cloth_origin = wp.vec3(0.6, 0.0, 0.102)

        # Cloth parameters: user's current material baseline
        self.density = 1.0
        self.tri_ke = 8.0e2 * 4
        self.tri_ka = 8.0e2 * 4
        self.tri_kd = 1.2e-1
        self.edge_ke = 6.0
        self.edge_kd = 6.0e-1
        self.particle_radius = 0.002

        # Grasp primitive parameters
        self.patch_side = 1                  # fixed single-point corner topology patch
        self.diag_patch_side = 3             # larger local region for diagnostics
        self.max_patch_step = 0.0008
        self.max_patch_step_lift = 0.00015   # slow-lift branch diagnostic
        self.alpha_base = 0.25
        self.grasp_rate = 0.06
        self.release_rate = 0.08

        builder = newton.ModelBuilder(gravity=-9.81)

        vertices = []
        indices = []
        for i in range(self.cloth_res):
            for j in range(self.cloth_res):
                x = (i / (self.cloth_res - 1) - 0.5) * self.cloth_size
                y = (j / (self.cloth_res - 1) - 0.5) * self.cloth_size
                vertices.append(wp.vec3(x, y, 0.0))
        for i in range(self.cloth_res - 1):
            for j in range(self.cloth_res - 1):
                idx = i * self.cloth_res + j
                indices += [idx, idx + 1, idx + self.cloth_res + 1]
                indices += [idx, idx + self.cloth_res + 1, idx + self.cloth_res]

        self.indices_flat = indices

        builder.add_cloth_mesh(
            vertices=vertices,
            indices=indices,
            pos=self.cloth_origin,
            rot=wp.quat_identity(),
            vel=wp.vec3(0.0, 0.0, 0.0),
            scale=1.0,
            density=self.density,
            tri_ke=self.tri_ke,
            tri_ka=self.tri_ka,
            tri_kd=self.tri_kd,
            edge_ke=self.edge_ke,
            edge_kd=self.edge_kd,
            particle_radius=self.particle_radius,
        )

        table_cfg = newton.ModelBuilder.ShapeConfig()
        table_cfg.density = 0.0
        table_cfg.ke = 1.5e4
        table_cfg.kd = 3e-1
        table_cfg.mu = 0.7
        body_table = builder.add_body(xform=wp.transform(self.table_center, wp.quat_identity()), label="table")
        builder.add_shape_box(body_table, hx=self.table_half[0], hy=self.table_half[1], hz=self.table_half[2], cfg=table_cfg)
        builder.add_ground_plane(cfg=table_cfg)

        # Visual/kinematic grasp sphere
        self.sphere_radius = 0.006
        body_sphere = builder.add_body(
            xform=wp.transform(wp.vec3(0.6, 0.0, 0.16), wp.quat_identity()),
            label="grasp_marker",
        )
        sphere_cfg = newton.ModelBuilder.ShapeConfig()
        sphere_cfg.density = 0.0
        sphere_cfg.ke = 1.0e5
        sphere_cfg.kd = 1.0e-4
        sphere_cfg.mu = 0.0
        builder.add_shape_sphere(body_sphere, radius=self.sphere_radius, cfg=sphere_cfg)
        self.sphere_body_index = 1  # table is 0, sphere is 1

        builder.color(include_bending=True)
        self.model = builder.finalize()
        if hasattr(self.model, "edge_rest_angle"):
            self.model.edge_rest_angle.zero_()

        self.model.shape_contact_pair_count = 0
        self.model.particle_friction = 0.8
        self.model.soft_contact_ke = 1.5e4
        self.model.soft_contact_kd = 3e-1
        self.model.soft_contact_mu = 0.7

        self.state_0 = self.model.state()
        self.state_1 = self.model.state()
        self.control = self.model.control()

        self.collision_pipeline = newton.CollisionPipeline(self.model, soft_contact_margin=0.0025)
        self.contacts = self.collision_pipeline.contacts()

        self.solver = SolverVBD(
            self.model,
            iterations=32,
            integrate_with_external_rigid_solver=True,
            particle_enable_self_contact=True,
            particle_self_contact_radius=0.0015,
            particle_self_contact_margin=0.0020,
            particle_topological_contact_filter_threshold=1,
            particle_rest_shape_contact_exclusion_radius=0.5,
            rigid_contact_k_start=2e4,
        )

        self.pristine_particle_q = self.state_0.particle_q.numpy().copy()
        self.pristine_inv_mass = self.model.particle_inv_mass.numpy().copy()

        # Edge / triangle diagnostics
        self.edge_pairs = build_unique_edges_from_triangles(indices)
        self.tris = build_triangles(indices)
        rest = self.pristine_particle_q
        self.edge_rest_lengths = np.linalg.norm(
            rest[self.edge_pairs[:, 0]] - rest[self.edge_pairs[:, 1]], axis=1
        )
        self.tri_rest_area, self.tri_rest_normal_unit = triangle_rest_geometry(rest, self.tris)

        self.viewer.set_model(self.model)
        if hasattr(self.viewer, "set_camera"):
            self.viewer.set_camera(wp.vec3(1.2, -0.8, 0.6), pitch=-20.0, yaw=40.0)

        self.reset_episode()

    def reset_episode(self):
        self.state_0.particle_q = wp.array(self.pristine_particle_q.copy(), dtype=wp.vec3, device=self.state_0.particle_q.device)
        self.state_0.particle_qd.zero_()
        self.state_1.particle_q = self.state_0.particle_q
        self.state_1.particle_qd.zero_()
        self.model.particle_inv_mass = wp.array(self.pristine_inv_mass.copy(), dtype=wp.float32, device=self.model.particle_inv_mass.device)
        self.solver.rebuild_bvh(self.state_0)

        particles = self.state_0.particle_q.numpy()
        cloth_center_xy = particles[:, :2].mean(axis=0)
        guess_xy = np.array([0.65, 0.05], dtype=np.float32)
        dists_corner = np.linalg.norm(particles[:, :2] - guess_xy[None, :], axis=1)
        self.corner_idx = int(np.argmin(dists_corner))
        true_corner = particles[self.corner_idx, :2]
        diag_vec = cloth_center_xy - true_corner
        diag_len = np.linalg.norm(diag_vec)
        diag_dir = diag_vec / max(diag_len, 1e-8)
        bite_dist = 0.005
        self.grasp_xy = true_corner + diag_dir * bite_dist
        self.target_xy = np.array([0.55, -0.05], dtype=np.float32)

        # Fixed topology grasp patch and a larger local diagnostic region around the corner
        self.topo_grasp_patch = build_corner_topology_patch(self.corner_idx, self.cloth_res, side=self.patch_side)
        self.topo_diag_patch = build_corner_topology_patch(self.corner_idx, self.cloth_res, side=self.diag_patch_side)
        self.target_eval_patch = select_target_patch_spatial(particles, self.target_xy, radius_xy=0.012, n_pts=6)

        diag_set = set(int(x) for x in self.topo_diag_patch.tolist())
        mask_edge = np.array([(int(i) in diag_set) or (int(j) in diag_set) for i, j in self.edge_pairs], dtype=bool)
        self.local_edge_pairs = self.edge_pairs[mask_edge]
        self.local_edge_rest_lengths = self.edge_rest_lengths[mask_edge]

        mask_tri = np.array([any(int(v) in diag_set for v in tri) for tri in self.tris], dtype=bool)
        self.local_tris = self.tris[mask_tri]
        self.local_tri_rest_area = self.tri_rest_area[mask_tri]
        self.local_tri_rest_normal_unit = self.tri_rest_normal_unit[mask_tri]

        self.max_global_stretch_seen = 1.0
        self.max_local_stretch_seen = 1.0
        self.min_global_area_seen = 1.0
        self.min_local_area_seen = 1.0
        self.max_global_flip_seen = 0
        self.max_local_flip_seen = 0

        self.phase = 0
        self.phase_timer = 0.0
        self.grasp_strength = 0.0
        self.is_grasped = False
        self.patch_idx = np.array([], dtype=np.int32)
        self.patch_local_offsets = np.zeros((0, 3), dtype=np.float32)
        self.patch_center = np.array([self.grasp_xy[0], self.grasp_xy[1], 0.16], dtype=np.float32)
        self.patch_target = self.patch_center.copy()
        self.release_gap = None
        self.finished = False

        body_q = self.state_0.body_q.numpy()
        body_q[self.sphere_body_index][0] = self.patch_center[0]
        body_q[self.sphere_body_index][1] = self.patch_center[1]
        body_q[self.sphere_body_index][2] = self.patch_center[2]
        self.state_0.body_q = wp.array(body_q, dtype=wp.transform)
        self.state_1.body_q = self.state_0.body_q

        print("\n[BRANCH] Reset cloth-only diagonal fold experiment")
        print(f"[BRANCH] grasp_xy=({self.grasp_xy[0]:.4f}, {self.grasp_xy[1]:.4f}), target_xy=({self.target_xy[0]:.4f}, {self.target_xy[1]:.4f})")
        print(f"[BRANCH] fixed topo patch = {self.topo_grasp_patch.tolist()}")

    def _update_kinematic_sphere(self):
        body_q = self.state_0.body_q.numpy()
        body_q[self.sphere_body_index][0] = self.patch_center[0]
        body_q[self.sphere_body_index][1] = self.patch_center[1]
        body_q[self.sphere_body_index][2] = self.patch_center[2]
        self.state_0.body_q = wp.array(body_q, dtype=wp.transform)
        self.state_1.body_q = self.state_0.body_q

    def _establish_grasp(self):
        particles = self.state_0.particle_q.numpy()
        patch_idx = self.topo_grasp_patch.copy()
        if len(patch_idx) == 0:
            return False
        self.patch_idx = patch_idx
        patch_pts = particles[patch_idx]
        patch_center0 = patch_pts.mean(axis=0)
        self.patch_center = patch_center0.copy()
        self.patch_target = patch_center0.copy()
        self.patch_local_offsets = patch_pts - patch_center0
        self.is_grasped = True
        self.grasp_strength = 0.0
        print(f"[BRANCH] grasp established with fixed topo patch of {len(self.patch_idx)} points at z={patch_center0[2]:.4f}")
        return True

    def _apply_patch_control(self):
        # Always move the grasp marker / target center, even before grasp is established.
        delta = self.patch_target - self.patch_center
        dn = np.linalg.norm(delta)

        step_limit = self.max_patch_step_lift if self.phase == 3 else self.max_patch_step
        if dn > step_limit:
            delta = delta / dn * step_limit
        self.patch_center = self.patch_center + delta

        # Only start driving cloth particles after grasp has been established.
        if not self.is_grasped or len(self.patch_idx) == 0:
            return

        particles = self.state_0.particle_q.numpy()
        vels = self.state_0.particle_qd.numpy()

        self.grasp_strength = min(1.0, self.grasp_strength + self.grasp_rate)
        alpha = self.alpha_base * self.grasp_strength
        vel_damp = 1.0 - 0.7 * self.grasp_strength

        target_pts = self.patch_center[None, :] + self.patch_local_offsets
        for i, idx in enumerate(self.patch_idx):
            particles[idx] = (1.0 - alpha) * particles[idx] + alpha * target_pts[i]
            vels[idx] *= vel_damp

        self.state_0.particle_q = wp.array(particles, dtype=wp.vec3, device=self.state_0.particle_q.device)
        self.state_0.particle_qd = wp.array(vels, dtype=wp.vec3, device=self.state_0.particle_qd.device)

        inv = self.pristine_inv_mass.copy()
        inv[self.patch_idx] = self.pristine_inv_mass[self.patch_idx] * (1.0 - 0.95 * self.grasp_strength)
        self.model.particle_inv_mass = wp.array(inv, dtype=wp.float32, device=self.model.particle_inv_mass.device)

    def _apply_release(self):
        if not self.is_grasped:
            return
        self.grasp_strength = max(0.0, self.grasp_strength - self.release_rate)
        inv = self.pristine_inv_mass.copy()
        if len(self.patch_idx) > 0:
            inv[self.patch_idx] = self.pristine_inv_mass[self.patch_idx] * (1.0 - 0.95 * self.grasp_strength)
        self.model.particle_inv_mass = wp.array(inv, dtype=wp.float32, device=self.model.particle_inv_mass.device)
        if self.grasp_strength <= 1e-4:
            self.is_grasped = False
            self.patch_idx = np.array([], dtype=np.int32)
            self.patch_local_offsets = np.zeros((0, 3), dtype=np.float32)
            self.model.particle_inv_mass = wp.array(self.pristine_inv_mass.copy(), dtype=wp.float32, device=self.model.particle_inv_mass.device)

    def _phase_logic(self):
        self.phase_timer += self.sim_dt

        if self.phase == 0:
            self.patch_target = np.array([self.grasp_xy[0], self.grasp_xy[1], 0.16], dtype=np.float32)
            if self.phase_timer > 0.5:
                self.phase = 1
                self.phase_timer = 0.0
                print("[BRANCH] phase 1: descend")

        elif self.phase == 1:
            self.patch_target = np.array([self.grasp_xy[0], self.grasp_xy[1], 0.108], dtype=np.float32)
            if self.patch_center[2] <= 0.112:
                ok = self._establish_grasp()
                self.phase = 2 if ok else 8
                self.phase_timer = 0.0
                print("[BRANCH] phase 2: engage" if ok else "[BRANCH] grasp failed")

        elif self.phase == 2:
            self.patch_target = np.array([self.grasp_xy[0], self.grasp_xy[1], 0.108], dtype=np.float32)
            if self.phase_timer > 0.20:
                self.phase = 3
                self.phase_timer = 0.0
                print("[BRANCH] phase 3: lift")

        elif self.phase == 3:
            self.patch_target = np.array([self.grasp_xy[0], self.grasp_xy[1], 0.205], dtype=np.float32)
            if self.patch_center[2] >= 0.200:
                self.phase = 4
                self.phase_timer = 0.0
                print("[BRANCH] phase 4: pause after lift (slow-lift diagnostic)")

        elif self.phase == 4:
            self.patch_target = np.array([self.grasp_xy[0], self.grasp_xy[1], 0.205], dtype=np.float32)
            if self.phase_timer > 0.50:
                self.finished = True

        elif self.phase == 5:
            self.patch_target = np.array([self.target_xy[0], self.target_xy[1], 0.205], dtype=np.float32)
            if np.linalg.norm(self.patch_center[:2] - self.target_xy) < 0.003:
                self.phase = 6
                self.phase_timer = 0.0
                print("[BRANCH] phase 6: pause over target")

        elif self.phase == 6:
            self.patch_target = np.array([self.target_xy[0], self.target_xy[1], 0.205], dtype=np.float32)
            if self.phase_timer > 0.20:
                particles = self.state_0.particle_q.numpy()
                if len(self.patch_idx) > 0:
                    grasp_xy = particles[self.patch_idx][:, :2].mean(axis=0)
                    target_xy = particles[self.target_eval_patch][:, :2].mean(axis=0)
                    self.release_gap = float(np.linalg.norm(grasp_xy - target_xy))
                self.phase = 7
                self.phase_timer = 0.0
                print("[BRANCH] phase 7: release")

        elif self.phase == 7:
            self.patch_target = np.array([self.target_xy[0], self.target_xy[1], 0.205], dtype=np.float32)
            self._apply_release()
            if not self.is_grasped:
                self.phase = 8
                self.phase_timer = 0.0
                print("[BRANCH] phase 8: settle")

        elif self.phase == 8:
            self.patch_target = np.array([self.target_xy[0], self.target_xy[1], 0.25], dtype=np.float32)
            if self.phase_timer > 0.8:
                self.finished = True

    def simulate(self):
        for _ in range(self.sim_substeps):
            self.state_0.clear_forces()
            self.viewer.apply_forces(self.state_0)

            self._phase_logic()
            self._apply_patch_control()
            self._update_kinematic_sphere()

            self.collision_pipeline.collide(self.state_0, self.contacts)
            self.solver.step(self.state_0, self.state_1, self.control, self.contacts, self.sim_dt)
            self.state_0, self.state_1 = self.state_1, self.state_0

            p = self.state_0.particle_q.numpy()
            gmax, _, _, _ = stretch_stats(p, self.edge_pairs, self.edge_rest_lengths)
            lmax, _, _, _ = stretch_stats(p, self.local_edge_pairs, self.local_edge_rest_lengths)
            gmin_area, _, _, gflips, _ = triangle_signed_area_stats(
                p, self.tris, self.tri_rest_area, self.tri_rest_normal_unit
            )
            lmin_area, _, _, lflips, _ = triangle_signed_area_stats(
                p, self.local_tris, self.local_tri_rest_area, self.local_tri_rest_normal_unit
            )

            self.max_global_stretch_seen = max(self.max_global_stretch_seen, gmax)
            self.max_local_stretch_seen = max(self.max_local_stretch_seen, lmax)
            self.min_global_area_seen = min(self.min_global_area_seen, gmin_area)
            self.min_local_area_seen = min(self.min_local_area_seen, lmin_area)
            self.max_global_flip_seen = max(self.max_global_flip_seen, gflips)
            self.max_local_flip_seen = max(self.max_local_flip_seen, lflips)

        self.sim_time += self.frame_dt

    def step(self):
        if not self.finished:
            self.simulate()
        else:
            self.report_and_reset()

    def report_and_reset(self):
        particles = self.state_0.particle_q.numpy()
        min_pos = particles.min(axis=0)
        max_pos = particles.max(axis=0)
        bbox_size = np.linalg.norm(max_pos - min_pos)
        max_vel = float(np.max(np.linalg.norm(self.state_0.particle_qd.numpy(), axis=1)))

        if len(self.patch_idx) > 0:
            grasp_xy = particles[self.patch_idx][:, :2].mean(axis=0)
        else:
            grasp_xy = particles[self.corner_idx, :2]
        target_xy = particles[self.target_eval_patch][:, :2].mean(axis=0)
        final_gap = float(np.linalg.norm(grasp_xy - target_xy))
        rebound = final_gap / max(self.release_gap or final_gap, 1e-6)

        gmax, gp95, gmean, _ = stretch_stats(particles, self.edge_pairs, self.edge_rest_lengths)
        lmax, lp95, lmean, _ = stretch_stats(particles, self.local_edge_pairs, self.local_edge_rest_lengths)

        garea_min, garea_p05, garea_mean, gflips, _ = triangle_signed_area_stats(
            particles, self.tris, self.tri_rest_area, self.tri_rest_normal_unit
        )
        larea_min, larea_p05, larea_mean, lflips, _ = triangle_signed_area_stats(
            particles, self.local_tris, self.local_tri_rest_area, self.local_tri_rest_normal_unit
        )

        print("\n[BRANCH][RESULT]")
        print(f"  final_gap             = {final_gap:.4f} m")
        print(f"  release_gap           = {(self.release_gap or float('nan')):.4f} m")
        print(f"  rebound_ratio         = {rebound:.3f}")
        print(f"  max_particle_v        = {max_vel:.4f} m/s")
        print(f"  bbox_size             = {bbox_size:.4f} m")
        print(f"  z_min                 = {min_pos[2]:.4f} m")
        print(f"  global_stretch_max    = {gmax:.3f}x")
        print(f"  global_stretch_p95    = {gp95:.3f}x")
        print(f"  global_stretch_mean   = {gmean:.3f}x")
        print(f"  local_stretch_max     = {lmax:.3f}x")
        print(f"  local_stretch_p95     = {lp95:.3f}x")
        print(f"  local_stretch_mean    = {lmean:.3f}x")
        print(f"  peak_global_seen      = {self.max_global_stretch_seen:.3f}x")
        print(f"  peak_local_seen       = {self.max_local_stretch_seen:.3f}x")
        print(f"  global_area_min       = {garea_min:.3f}x")
        print(f"  global_area_p05       = {garea_p05:.3f}x")
        print(f"  global_area_mean      = {garea_mean:.3f}x")
        print(f"  global_flip_count     = {gflips}")
        print(f"  local_area_min        = {larea_min:.3f}x")
        print(f"  local_area_p05        = {larea_p05:.3f}x")
        print(f"  local_area_mean       = {larea_mean:.3f}x")
        print(f"  local_flip_count      = {lflips}")
        print(f"  peak_min_global_area  = {self.min_global_area_seen:.3f}x")
        print(f"  peak_min_local_area   = {self.min_local_area_seen:.3f}x")
        print(f"  peak_global_flips     = {self.max_global_flip_seen}")
        print(f"  peak_local_flips      = {self.max_local_flip_seen}")

        self.reset_episode()

    def render(self):
        self.viewer.begin_frame(self.sim_time)
        self.viewer.log_state(self.state_0)
        self.viewer.end_frame()


if __name__ == "__main__":
    parser = newton.examples.create_parser()
    viewer, args = newton.examples.init(parser)
    example = Example(viewer, args)
    newton.examples.run(example, args)

[jumyungc]

Thank you for reporting this!

I'm guessing the stability issue is the most critical one:

if I directly drive cloth particles or a small patch, the cloth often becomes numerically unstable during lifting, with severe local distortion and unrealistic deformation;

I'd recommend tuning all kd values down to something like 1.0e-4 and increasing from there as needed.
In the current VBD damping convention, the effective damping is actually ke * kd (not kd alone), so the values are much larger than they appear. Setting damping too high causes the instability you observed.
Higher damping will stabilize behavior, but there's a threshold beyond which it becomes unstable. We are currently working on making this much easier to use!

I'd also highly recommend using CUDA graph capture, which you are not currently using. It will improve performance significantly.

One minor note: you may not need to set rigid_contact_k_start, or if you do, just use the same value as the table contact stiffness to keep tuning simple. We are working on improving this as well, so please stay tuned.

Please let us know if you encounter any other issues.

[N3rv3rM1nd]

Thank you so much for the detailed response and valuable suggestions! I really appreciate you clarifying the ke * kd damping convention and providing these practical tips to address the stability issue.

I am currently running new experiments following your guidance exactly. I am tuning the kd values starting from 1.0e-4, looking into implementing CUDA graph capture to improve performance, and simplifying the rigid_contact_k_start setting as you recommended.

I will let you know how these adjustments work out or if I encounter any further issues. Thank you again for your time and your ongoing efforts to improve the system!

Best regards,