new line search for strain limiting implemented

Author: ryichando

README.md

@@ -28,6 +28,7 @@ Published in [ACM Transactions on Graphics (TOG)](https://dl.acm.org/doi/abs/10.
 - [⚡️ Requirements](#️-requirements)
 - [📝 Change History](#-change-history)
 - [🐍 How To Use](#-how-to-use)
+- [🔍 Obtaining Logs](#-obtaining-logs)
 - [🖼️ Catalogue](#️-catalogue)
 - [🚀 GitHub Actions](#-github-actions)
 - [💨 Getting Started](#-getting-started)
@@ -64,6 +65,7 @@ Published in [ACM Transactions on Graphics (TOG)](https://dl.acm.org/doi/abs/10.
 
 ## 📝 Change History
 
+- (2024.12.27) Line search for strain limiting is improved [[Markdown]](./articles/bug.md#new-strain-limiting-line-search).
 - (2024.12.23) Added [[Bug Fixes and Updates]](./articles/bug.md)
 - (2024.12.21) Added a [house of cards example](./examples/cards.ipynb) [[Video]](https://drive.google.com/file/d/1PMdDnlyCsjinbvICKph_0UcXUfUvvUmZ/view)
 - (2024.12.18) Added a [frictional contact example](./examples/friction.ipynb): armadillo sliding on the slope [[Video]](https://drive.google.com/file/d/12WGdfDTFIwCT0UFGEZzfmQreM6WSSHet/view)
@@ -164,13 +166,13 @@ session.export.animation(path)
 ```
 <img src="./asset/image/drape.jpg" alt="drape">
 
-### 🔍 Obtaining Logs
+## 🔍 Obtaining Logs
 
 Logs for the simulation can also be queried through the Python APIs. Here's an example of how to get the list of recorded logs, fetch them, and compute the average.
 
 ```python
-# get list of logs files list[str]
-logs = session.get.log()
+# get the list of log names as list[str]
+logs = session.get.logfiles()
 assert 'per_video_frame' in logs
 assert 'advance.newton_steps' in logs
 
@@ -187,23 +189,48 @@ newton_steps = session.get.numbers('advance.newton_steps')
 print('avg:', sum([n for _,n in newton_steps])/len(newton_steps))
 ```
 
-Here are the representative ones.
+Below are some representatives.
 `vid_time` refers to the video time in seconds and is recorded as `float`.
-`msec` refers to the consumed simulation time recorded as `int`.
+`ms` refers to the consumed simulation time in milliseconds recorded as `int`.
 `vid_frame` is the video frame count recorede as `int`.
 
 | **Log Name** | **Description** | **Format**
 |---------------|----------------|------------
-| **per_video_frame** | Time per video frame | list[(vid_frame,msec)] |
-| **advance.matrix_assembly** | Matrix assembly time | list[(vid_time,msec)] |
-| **advance.linsolve** | Linear system solve time | list[(vid_time,msec)] |
-| **advance.line_search** | Line search time | list[(vid_time,msec)] |
-| **advance** | Time per step | list[(vid_time,msec)] |
+| **per_video_frame** | Time per video frame | list[(vid_frame,ms)] |
+| **advance.matrix_assembly** | Matrix assembly time | list[(vid_time,ms)] |
+| **advance.linsolve** | Linear system solve time | list[(vid_time,ms)] |
+| **advance.line_search** | Line search time | list[(vid_time,ms)] |
+| **advance** | Time per step | list[(vid_time,ms)] |
 | **advance.newton_steps** | Newton iterations per step | list[(vid_time,count)] |
 | **advance.num_contact** | Contact count | list[(vid_time,count)] |
 | **advance.max_sigma** | Max stretch | list(vid_time,strech) |
 
-Note that some entries have multiple records at the same video time ⏱️. This occurs because the same operation is executed multiple times 🔄 within a single step during the inner Newton's iteration 🧮. For example, the linear system solve is performed at each Newton's step, so if multiple Newton's steps are 🔁 executed, multiple linear system solve times may appear in the record at the same 📊 video time.
+Note that some entries have multiple records at the same video time ⏱️. This occurs because the same operation is executed multiple times 🔄 within a single step during the inner Newton's iterations 🧮. For example, the linear system solve is performed at each Newton's step, so if multiple Newton's steps are 🔁 executed, multiple linear system solve times may appear in the record at the same 📊 video time.
+
+If you would like to retrieve the raw log stream, you can do so by
+
+```python
+# Last 8 lines. Omit for everything.
+for line in session.get.log(n_lines=8):
+    print(line)
+```
+
+This will output something like:
+
+```
+* dt: 1.000e-03
+* max_sigma: 1.045e+00
+* avg_sigma: 1.030e+00
+------ newton step 1 ------
+   ====== contact_matrix_assembly ======
+   > dry_pass...0 msec
+   > rebuild...7 msec
+   > fillin_pass...0 msec
+```
+
+If you would like to read `stdout` and `stderr`, you can do so using `session.get.stdout()` and `session.get.stderr()` (if it exists). They return `list[str]`.
+
+All the log files 📂 are available ✅ and can be fetched ⬇️ during the simulation 🧑‍💻.
 
 ## 🖼️ Catalogue
 
@@ -568,22 +595,22 @@ The author also extends thanks to the teams in the IP department for permitting
 
 ```
 @article{Ando2024CB,
-author = {Ando, Ryoichi},
-title = {A Cubic Barrier with Elasticity-Inclusive Dynamic Stiffness},
-year = {2024},
-issue_date = {December 2024},
-publisher = {Association for Computing Machinery},
-address = {New York, NY, USA},
-volume = {43},
-number = {6},
-issn = {0730-0301},
-url = {https://doi.org/10.1145/3687908},
-doi = {10.1145/3687908},
-journal = {ACM Trans. Graph.},
-month = nov,
-articleno = {224},
-numpages = {13},
-keywords = {collision, contact}
+    author = {Ando, Ryoichi},
+    title = {A Cubic Barrier with Elasticity-Inclusive Dynamic Stiffness},
+    year = {2024},
+    issue_date = {December 2024},
+    publisher = {Association for Computing Machinery},
+    address = {New York, NY, USA},
+    volume = {43},
+    number = {6},
+    issn = {0730-0301},
+    url = {https://doi.org/10.1145/3687908},
+    doi = {10.1145/3687908},
+    journal = {ACM Trans. Graph.},
+    month = nov,
+    articleno = {224},
+    numpages = {13},
+    keywords = {collision, contact}
 }
 ```
 

articles/bug.md

@@ -1,12 +1,53 @@
 ## Bug Fixes 🐞 and Updates 🔄
 
+### New Strain Limiting Line Search
+
+On 2024-Dec-27, we made minor improvements to the line search algorithm for strain limiting.
+The new algorithm is laid out below.
+Compared to the one presented in the supplementary PDF, this new algorithm searches for `t` (time of impact) with exact numerical precision; no parameters, such as the `eps` tolerance for checking convergence or the maximum loop count, are necessary.
+This new approach converges quickly and exits the loop fast.
+
+```c++
+/* Apache v2.0 License */
+float line_search_strain_limiting(
+    const Mat3x2f &F0, // Deformation gradient at t = 0.0
+    const Mat3x2f &F1, // Deformation gradient at t = 1.0
+    float t,           // Maximal time of impact, such as t = 1.0
+    float max_sigma )  // Maximal sigma, such as 1.01
+{
+    const Mat3x2f dF = F1 - F0;
+    if (svd3x2(F0 + t * dF).S.max() >= max_sigma) {
+        float upper_t = t;
+        float lower_t = 0.0f;
+        float window = upper_t - lower_t;
+        while (true) {
+            t = 0.5f * (upper_t + lower_t);
+            float diff = svd3x2(F0 + t * dF).S.maxCoeff() - max_sigma;
+            if (diff < 0.0f) {
+                lower_t = t;
+            } else {
+                upper_t = t;
+            }
+            float new_window = upper_t - lower_t;
+            if (new_window == window) {
+                break;
+            } else {
+                window = new_window;
+            }
+        }
+        t = lower_t;
+    }
+    return t;
+}
+```
+
 ### BVH Construction
 
 In the supplementary PDF, we mentioned that the BVH is reconstructed every 10 video frames; however, this is not what is implemented in the public code.
 
 In this code, the BVH is continuously updated on the CPU side in the background without blocking the simulation. At the beginning of each step, we check if the BVH construction is finished, and if it is, we update the GPU buffer.
 
-### Hang Example (2024-Dec-22)
+### Hang Example
 
 On 2024-Dec-22, we encountered a situation where the PCG solver in our "hang" example `hang.ipynb` failed to converge, resulting in a simulation failure.
 

examples/frontend/_session_.py

@@ -266,14 +266,40 @@ class SessionGet:
     def __init__(self, session: "Session"):
         self._session = session
 
-    def log(self) -> list[str]:
+    def logfiles(self) -> list[str]:
         path = os.path.join(self._session.info.path, "output", "data")
         result = []
         for file in os.listdir(path):
             if file.endswith(".out"):
                 result.append(file.replace(".out", ""))
         return result
 
+    def _tail_file(self, path: str, n_lines: Optional[int] = None) -> list[str]:
+        if os.path.exists(path):
+            with open(path, "r") as f:
+                lines = f.readlines()
+                lines = [line.rstrip("\n") for line in lines]
+                if n_lines is not None:
+                    return lines[-n_lines:]
+                else:
+                    return lines
+        return []
+
+    def log(self, n_lines: Optional[int] = None) -> list[str]:
+        return self._tail_file(
+            os.path.join(self._session.info.path, "output", "cudasim_log.txt"), n_lines
+        )
+
+    def stdout(self, n_lines: Optional[int] = None) -> list[str]:
+        return self._tail_file(
+            os.path.join(self._session.info.path, "stdout.log"), n_lines
+        )
+
+    def stderr(self, n_lines: Optional[int] = None) -> list[str]:
+        return self._tail_file(
+            os.path.join(self._session.info.path, "error.log"), n_lines
+        )
+
     def numbers(self, name: str):
         def float_or_int(var):
             var = float(var)
@@ -427,13 +453,10 @@ def init(self, scene: FixedScene) -> "Session":
 
     def finished(self) -> bool:
         finished_path = os.path.join(self.output.path, "finished.txt")
-        err_path = os.path.join(self.info.path, "error.log")
-        if os.path.exists(err_path):
-            file = open(err_path, "r")
-            lines = file.readlines()
-            if len(lines) > 0:
-                for line in lines:
-                    print(line)
+        error = self.get.stderr()
+        if len(error) > 0:
+            for line in error:
+                print(line)
         return os.path.exists(finished_path)
 
     def start(self, param: Param, force: bool = True, blocking=False) -> "Session":

src/backend.rs

@@ -133,7 +133,7 @@ impl Backend {
                 let face = builder::build_bvh(&vertex, Some(&face));
                 let edge = builder::build_bvh(&vertex, Some(&edge));
                 let vertex = builder::build_bvh::<1>(&vertex, None);
-                result_sender.send(BvhSet { face, edge, vertex }).unwrap();
+                let _ = result_sender.send(BvhSet { face, edge, vertex });
             }
         });
 

src/cpp/strainlimiting/strainlimiting.cu

@@ -109,26 +109,27 @@ float line_search(const DataSet &data, const Kinematic &kinematic,
                 utility::compute_deformation_grad(x1, inv_rest2x2[i]);
             const Mat3x2f dF = F1 - F0;
             if (utility::svd3x2(F0 + t * dF).S.maxCoeff() >= max_sigma) {
-                float tol = param.ccd_reduction *
-                            (max_sigma - utility::svd3x2(F0).S.maxCoeff());
-                float target = max_sigma - 2.0f * tol;
                 float upper_t = t;
                 float lower_t = 0.0f;
-                for (unsigned k = 0; k < param.binary_search_max_iter; ++k) {
+                float window = upper_t - lower_t;
+                while (true) {
                     t = 0.5f * (upper_t + lower_t);
                     Svd3x2 svd = utility::svd3x2(F0 + t * dF);
-                    float diff = svd.S.maxCoeff() - target;
-                    if (fabs(diff) < tol) {
-                        break;
-                    } else if (diff < 0.0f) {
+                    float diff = svd.S.maxCoeff() - max_sigma;
+                    if (diff < 0.0f) {
                         lower_t = t;
                     } else {
                         upper_t = t;
                     }
+                    float new_window = upper_t - lower_t;
+                    if (new_window == window) {
+                        break;
+                    } else {
+                        window = new_window;
+                    }
                 }
-                while (utility::svd3x2(F0 + t * dF).S.maxCoeff() >= target) {
-                    t *= param.dt_decrease_factor;
-                }
+                t = lower_t;
+                assert(t > std::numeric_limits<float>::epsilon());
             }
         }
         toi[i] = t;