self-contact find boundary

Author: liminchen

self_contact/readme.md

@@ -1,6 +1,6 @@
 # Inversion-free Hyperelastic Solids Simulation
 
-Two squares falling onto the ground under gravity, contacting with each other, is simulated with an inversion-free hyperelastic potential and IPC with implicit Euler time integration.
+Two squares falling onto the ground under gravity, contacting with each other, and then compressed by a ceiling is simulated with an inversion-free hyperelastic potential and IPC with implicit Euler time integration.
 Each time step is solved by minimizing the Incremental Potential with the projected Newton method.
 
 ## Dependencies

self_contact/simulator.py

@@ -1,4 +1,4 @@
-# Mass-Spring Solids Simulation
+# FEM Solids Simulation
 
 import numpy as np  # numpy for linear algebra
 import pygame       # pygame for visualization
@@ -26,6 +26,7 @@
 [x, e] = square_mesh.generate(side_len, n_seg)       # node positions and triangle node indices of the top square
 e = np.append(e, np.array(e) + [len(x)] * 3, axis=0) # add triangle node indices of the bottom square
 x = np.append(x, x + [side_len * 0.1, -side_len * 1.1], axis=0) # add node positions of the bottom square
+[bp, be] = square_mesh.find_boundary(e)             # find boundary points and edges for self-contact
 x = np.append(x, [[0.0, side_len * 0.6]], axis=0)   # ceil origin (with normal [0.0, -1.0])
 v = np.array([[0.0, 0.0]] * len(x))                 # velocity
 m = [rho * side_len * side_len / ((n_seg + 1) * (n_seg + 1))] * len(x)  # calculate node mass evenly

self_contact/square_mesh.py

@@ -20,4 +20,23 @@ def generate(side_length, n_seg):
                 e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1])
                 e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1, i * (n_seg + 1) + j + 1])
 
-    return [x, e]
+    return [x, e]
+
+def find_boundary(e):
+    # index all half-edges for fast query
+    edge_set = set()
+    for i in range(0, len(e)):
+        for j in range(0, 3):
+            edge_set.add((e[i][j], e[i][(j + 1) % 3]))
+
+    # find boundary points and edges
+    bp_set = set()
+    be = []
+    for eI in edge_set:
+        if (eI[1], eI[0]) not in edge_set:
+            # if the inverse edge of a half-edge does not exist,
+            # then it is a boundary edge
+            be.append([eI[0], eI[1]])
+            bp_set.add(eI[0])
+            bp_set.add(eI[1])
+    return [list(bp_set), be]