trapped example with squishy ball

examples/trapped.ipynb

@@ -7,17 +7,44 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "import os\n",
+    "\n",
+    "import numpy as np\n",
+    "import sdf\n",
     "from frontend import App\n",
     "\n",
     "app = App.create(\"trapped\")\n",
+    "filepath = \"/tmp/squishy.tmp.ply\"\n",
     "\n",
-    "V, F, T = app.mesh.preset(\"armadillo\").decimate(30000).tetrahedralize().normalize()\n",
-    "app.asset.add.tet(\"armadillo\", V, F, T)\n",
+    "# create squishy ball mesh if not exists\n",
+    "if not os.path.exists(filepath):\n",
+    "    V, F = app.mesh.icosphere(r=2, subdiv_count=2)\n",
+    "    func = sdf.sphere(1.1)\n",
+    "    for f in F:\n",
+    "        d = np.mean(V[f], axis=0)\n",
+    "        if d[0] > 0:\n",
+    "            func = func | sdf.capsule(-d, d, 0.05)\n",
+    "    func.save(filepath, step=0.03)\n",
     "\n",
-    "scene = app.scene.create(\"sphere-trap\")\n",
+    "V, F, T = (\n",
+    "    app.mesh.load_tri(filepath)\n",
+    "    .decimate(100000)\n",
+    "    .tetrahedralize()\n",
+    "    .normalize()\n",
+    "    .scale(0.97)\n",
+    ")\n",
+    "app.asset.add.tet(\"squishy\", V, F, T)\n",
     "\n",
-    "scene.add.invisible.sphere([0, 0, 0], 0.7).invert().radius(0.15, 3).radius(100, 4)\n",
-    "scene.add(\"armadillo\").jitter().rotate(180, \"y\")\n",
+    "scene = app.scene.create(\"sphere-trap\")\n",
+    "(\n",
+    "    scene.add.invisible.sphere([0, 0, 0], 0.98)\n",
+    "    .invert()\n",
+    "    .radius(0.4, 2)\n",
+    "    .radius(0.4, 3)\n",
+    "    .radius(10, 4)\n",
+    ")\n",
+    "scene.add(\"squishy\").at(0.5, 0, 0).jitter()\n",
+    "scene.add(\"squishy\").at(-0.5, 0, 0).jitter()\n",
     "\n",
     "fixed = scene.build().report()\n",
     "fixed.preview()"
@@ -31,7 +58,13 @@
    "outputs": [],
    "source": [
     "param = app.session.param()\n",
-    "param.set(\"gravity\", 0.0).set(\"dt\", 0.01)\n",
+    "(\n",
+    "    param.set(\"gravity\", 0.0)\n",
+    "    .set(\"friction\", 0.0)\n",
+    "    .set(\"csrmat-max-nnz\", 3000000)\n",
+    "    .set(\"dt\", 0.01)\n",
+    ")\n",
+    "param.dyn(\"playback\").time(2.99).hold().time(3).change(0.1)\n",
     "\n",
     "session = app.session.create(fixed)\n",
     "session.start(param).preview()\n",

frontend/_mesh_.py

@@ -198,8 +198,6 @@ def icosphere(self, r: float = 1, subdiv_count: int = 3) -> "TriMesh":
 
     def _from_o3d(self, o3d_mesh) -> "TriMesh":
         """Load a mesh from an Open3D mesh"""
-        if o3d_mesh.is_self_intersecting():
-            print("Warning: Mesh is self-intersecting")
         return TriMesh.create(
             np.asarray(o3d_mesh.vertices),
             np.asarray(o3d_mesh.triangles),
@@ -819,16 +817,15 @@ def recompute_hash(self) -> "TriMesh":
         """Recompute the hash of the mesh"""
         import hashlib
 
-        self.hash = hashlib.sha256(
-            np.concatenate(
-                [
-                    np.array(self[0].shape),
-                    self[0].ravel(),
-                    np.array(self[1].shape),
-                    self[1].ravel(),
-                ]
-            )
-        ).hexdigest()
+        obj = np.concatenate(
+            [
+                np.array(self[0].shape),
+                self[0].ravel(),
+                np.array(self[1].shape),
+                self[1].ravel(),
+            ]
+        )
+        self.hash = hashlib.sha256(obj.tobytes()).hexdigest()
         return self
 
     def set_cache_dir(self, cache_dir: str) -> "TriMesh":

frontend/_parse_.py

@@ -201,6 +201,8 @@ def clear_doc():
                                 default_value = (
                                     int(float_value)
                                     if float_value.is_integer()
+                                    and var_type != "f32"
+                                    and var_type != "f64"
                                     else float_value
                                 )
                             except ValueError:

frontend/_scene_.py

@@ -999,9 +999,7 @@ def preview(
 
             if len(self._rod):
                 if plot is None:
-                    plot = self._plot.create().curve(
-                        vert, self._rod, shading=shading
-                    )   
+                    plot = self._plot.create().curve(vert, self._rod, shading=shading)
                 else:
                     plot.add.edge(vert, self._rod)
 
@@ -1181,7 +1179,7 @@ def min(self, axis: str) -> float:
                 result = min(result, np.min(vert[:, _axis[axis]]))
         return result
 
-    def max(self, axis: int) -> float:
+    def max(self, axis: str) -> float:
         """Get the maximum value of the scene along a specific axis.
 
         Args:
@@ -1240,8 +1238,8 @@ def add_entry(
         for name, obj in self._object.items():
             if not obj.static and obj.get("T") is None:
                 map = tag[name]
-                vert, edge = obj.get("V"), obj.get("E")
-                if vert is not None and edge is not None:
+                edge = obj.get("E")
+                if edge is not None:
                     add_entry(
                         map,
                         edge,
@@ -1252,14 +1250,24 @@ def add_entry(
         for name, obj in self._object.items():
             if not obj.static and obj.get("T") is None:
                 map = tag[name]
-                vert, tri = obj.get("V"), obj.get("F")
-                if vert is not None and tri is not None:
+                tri = obj.get("F")
+                if tri is not None:
                     add_entry(
                         map,
                         tri,
                     )
         shell_vert_start, shell_vert_end = rod_vert_end, concat_count
 
+        for name, obj in self._object.items():
+            if not obj.static:
+                map = tag[name]
+                tri = obj.get("F")
+                if tri is not None:
+                    add_entry(
+                        map,
+                        tri,
+                    )
+
         pbar.update(1)
         for name, obj in self._object.items():
             if not obj.static:
@@ -1325,10 +1333,15 @@ def vec_map(map, elm):
         for name, obj in self._object.items():
             if not obj.static:
                 map = tag[name]
-                tet = obj.get("T")
                 tri = obj.get("F")
-                if tet is not None and tri is not None:
+                if tri is not None and obj.get("T") is not None:
                     concat_tri.extend(vec_map(map, tri))
+
+        for name, obj in self._object.items():
+            if not obj.static:
+                map = tag[name]
+                tet = obj.get("T")
+                if tet is not None:
                     concat_tet.extend(vec_map(map, tet))
 
         pbar.update(1)

frontend/_session_.py

@@ -136,10 +136,13 @@ def hold(self) -> "Param":
                 self.change(last_val)
             else:
                 val = self._param[self._key]["value"]
+                val_type = self._param[self._key]["type"]
                 if isinstance(val, float):
                     self.change(val)
                 else:
-                    raise ValueError("Key must be float")
+                    raise ValueError(
+                        f"Key must be float. {val} is given. type: {val_type}"
+                    )
         return self
 
     def export(self, path: str):
@@ -885,6 +888,7 @@ def start(self, param: Param, force: bool = True, blocking=False) -> "Session":
         while not os.path.exists(log_path) and not os.path.exists(err_path):
             time.sleep(1)
         if process.poll() is not None:
+            display_log(open(err_path, "r").readlines())
             raise ValueError("Solver failed to start")
         else:
             init_path = os.path.join(self.info.path, "output", "data", "initialize.out")
@@ -894,8 +898,7 @@ def start(self, param: Param, force: bool = True, blocking=False) -> "Session":
                     break
                 time.sleep(1)
             if not os.path.exists(init_path):
-                err_content = open(err_path, "r").readlines()
-                display_log(err_content)
+                display_log(open(err_path, "r").readlines())
                 raise ValueError("Solver failed to start")
         if blocking or not self._in_jupyter_notebook:
             print(f">>> Log path: {log_path}")
@@ -918,7 +921,7 @@ def start(self, param: Param, force: bool = True, blocking=False) -> "Session":
                 print("*** Solver FAILED ***")
             else:
                 print("*** Solver finished ***")
-            n_logs = 16
+            n_logs = 32
             log_lines = open(log_path, "r").readlines()
             print(">>> Log:")
             for line in log_lines[-n_logs:]:
@@ -992,21 +995,21 @@ def convert_integer(number) -> str:
                 elif number < 1000:
                     return str(number)
                 elif number < 1_000_000:
-                    return f"{number/1_000:.2f}k"
+                    return f"{number / 1_000:.2f}k"
                 elif number < 1_000_000_000:
-                    return f"{number/1_000_000:.2f}M"
+                    return f"{number / 1_000_000:.2f}M"
                 else:
-                    return f"{number/1_000_000_000:.2f}B"
+                    return f"{number / 1_000_000_000:.2f}B"
 
             def convert_time(time) -> str:
                 if time is None:
                     return "N/A"
                 elif time < 1_000:
                     return f"{int(time)}ms"
                 elif time < 60_000:
-                    return f"{time/1_000:.2f}s"
+                    return f"{time / 1_000:.2f}s"
                 else:
-                    return f"{time/60_000:.2f}m"
+                    return f"{time / 60_000:.2f}m"
 
             if live_update and is_running():
 
@@ -1243,3 +1246,6 @@ def display_log(lines: list[str]):
         text = "\n".join(lines)
         log_widget.value = CONSOLE_STYLE + f"<pre style='no-scroll'>{text}</pre>"
         display(log_widget)
+    else:
+        for line in lines:
+            print(line)

src/args.rs

@@ -28,6 +28,14 @@ pub struct Args {
     #[clap(long, default_value_t = 1e-3)]
     pub dt: f32,
 
+    // Name: Playback Speed
+    // Description:
+    // The speed at which the simulation is played back.
+    // 1.0 means unaltered playback, values greater than 1.0 speed up the playback,
+    // and values less than 1.0 slow down the playback.
+    #[clap(long, default_value_t = 1.0)]
+    pub playback: f32,
+
     // Name: Lower Bound of Newton Steps
     // Recommended Range: 0 to 32
     // Description:

src/builder.rs

@@ -340,6 +340,7 @@ pub fn make_param(args: &Args) -> data::ParamSet {
         constraint_ghat: args.constraint_ghat,
         prev_dt: dt,
         dt,
+        playback: args.playback,
         min_newton_steps: args.min_newton_steps,
         target_toi: args.target_toi,
         bend: args.bend,

src/cpp/data.hpp

@@ -218,6 +218,7 @@ struct ParamSet {
     float constraint_ghat;
     float prev_dt;
     float dt;
+    float playback;
     unsigned min_newton_steps;
     float target_toi;
     float bend;

src/cpp/main/main.cu

@@ -290,14 +290,20 @@ StepResult advance() {
         } DISPATCH_END;
     }
 
-    float dt = param->dt;
+    float dt = param->playback * param->dt;
 
     // Name: Step Size
     // Format: list[(vid_time,float)]
     // Description:
     // Target step size.
     logging.mark("dt", dt);
 
+    // Name: playback
+    // Format: list[(vid_time,float)]
+    // Description:
+    // Playback speed.
+    logging.mark("playback", param->playback);
+
     if (shell_face_count) {
         utility::compute_svd(data, data.vertex.curr, svd, prm);
         tmp_scalar.clear();
@@ -632,7 +638,7 @@ StepResult advance() {
                 logging.mark("final_dt", dt);
 
                 param->prev_dt = dt;
-                param->time += param->prev_dt;
+                param->time += param->prev_dt / param->playback;
 
                 dev_dataset.vertex.prev.copy(dev_dataset.vertex.curr);
                 dev_dataset.vertex.curr.copy(eval_x);

src/data.rs

@@ -292,6 +292,7 @@ pub struct ParamSet {
     pub constraint_ghat: f32,
     pub prev_dt: f32,
     pub dt: f32,
+    pub playback: f32,
     pub min_newton_steps: u32,
     pub target_toi: f32,
     pub bend: f32,

src/scene.rs

@@ -619,6 +619,7 @@ impl Scene {
                     match title.as_str() {
                         "gravity" => param.gravity = Vec3f::new(0.0, val, 0.0),
                         "dt" => param.dt = val,
+                        "playback" => param.playback = val,
                         "friction" => param.friction = val,
                         _ => (),
                     }