added anchors

Author: liminchen

7_self_contact/BarrierEnergy.py

@@ -24,6 +24,7 @@ def val(x, n, o, bp, be, contact_area):
         if d < dhat:
             s = d / dhat
             sum += contact_area[i] * dhat * kappa / 2 * (s - 1) * math.log(s)
+    # ANCHOR: value
     # self-contact
     dhat_sqr = dhat * dhat
     for xI in bp:
@@ -34,6 +35,7 @@ def val(x, n, o, bp, be, contact_area):
                     s = d_sqr / dhat_sqr
                     # since d_sqr is used, need to divide by 8 not 2 here for consistency to linear elasticity:
                     sum += 0.5 * contact_area[xI] * dhat * kappa / 8 * (s - 1) * math.log(s)
+    # ANCHOR_END: value
     return sum
 
 def grad(x, n, o, bp, be, contact_area):
@@ -53,6 +55,7 @@ def grad(x, n, o, bp, be, contact_area):
             local_grad = contact_area[i] * dhat * (kappa / 2 * (math.log(s) / dhat + (s - 1) / d)) * n
             g[i] += local_grad
             g[-1] -= local_grad
+    # ANCHOR: gradient
     # self-contact
     dhat_sqr = dhat * dhat
     for xI in bp:
@@ -66,6 +69,7 @@ def grad(x, n, o, bp, be, contact_area):
                     g[xI] += local_grad[0:2]
                     g[eI[0]] += local_grad[2:4]
                     g[eI[1]] += local_grad[4:6]
+    # ANCHOR_END: gradient
     return g
 
 def hess(x, n, o, bp, be, contact_area):
@@ -94,6 +98,7 @@ def hess(x, n, o, bp, be, contact_area):
                             IJV[0].append(index[nI] * 2 + r)
                             IJV[1].append(index[nJ] * 2 + c)
                             IJV[2] = np.append(IJV[2], ((-1) ** (nI != nJ)) * local_hess[r, c])
+    # ANCHOR: Hessian
     # self-contact
     dhat_sqr = dhat * dhat
     for xI in bp:
@@ -114,6 +119,7 @@ def hess(x, n, o, bp, be, contact_area):
                                     IJV[0].append(index[nI] * 2 + r)
                                     IJV[1].append(index[nJ] * 2 + c)
                                     IJV[2] = np.append(IJV[2], local_hess[nI * 2 + r, nJ * 2 + c])
+    # ANCHOR_END: Hessian
     return IJV
 
 def init_step_size(x, n, o, bp, be, p):
@@ -129,6 +135,7 @@ def init_step_size(x, n, o, bp, be, p):
         p_n = (p[i] - p[-1]).dot(n)
         if p_n < 0:
             alpha = min(alpha, 0.9 * n.dot(x[i] - x[-1]) / -p_n)
+    # ANCHOR: line_search_filtering
     # self-contact
     for xI in bp:
         for eI in be:
@@ -137,6 +144,7 @@ def init_step_size(x, n, o, bp, be, p):
                     toc = CCD.narrow_phase_CCD(x[xI], x[eI[0]], x[eI[1]], p[xI], p[eI[0]], p[eI[1]], alpha)
                     if alpha > toc:
                         alpha = toc
+    # ANCHOR_END: line_search_filtering
     return alpha
 
 def compute_mu_lambda(x, n, o, bp, be, contact_area, mu):

7_self_contact/distance/CCD.py

@@ -1,3 +1,4 @@
+# ANCHOR: broad_phase
 from copy import deepcopy
 import numpy as np
 import math
@@ -14,7 +15,9 @@ def bbox_overlap(p, e0, e1, dp, de0, de1, toc_upperbound):
         return False
     else:
         return True
+# ANCHOR_END: broad_phase
 
+# ANCHOR: accd
 # compute the first "time" of contact, or toc,
 # return the computed toc only if it is smaller than the previously computed toc_upperbound
 def narrow_phase_CCD(_p, _e0, _e1, _dp, _de0, _de1, toc_upperbound):
@@ -56,4 +59,5 @@ def narrow_phase_CCD(_p, _e0, _e1, _dp, _de0, _de1, toc_upperbound):
         if toc > toc_upperbound:
             return toc_upperbound
 
-    return toc
+    return toc
+# ANCHOR_END: accd

7_self_contact/distance/PointEdgeDistance.py

@@ -1,3 +1,4 @@
+# ANCHOR: PE_val_grad
 import numpy as np
 
 import distance.PointPointDistance as PP
@@ -24,6 +25,7 @@ def grad(p, e0, e1):
         return np.reshape([g_PP[0:2], np.array([0.0, 0.0]), g_PP[2:4]], (1, 6))[0]
     else:            # point(p)-line(e0e1) expression
         return PL.grad(p, e0, e1)
+# ANCHOR_END: PE_val_grad
 
 def hess(p, e0, e1):
     e = e1 - e0

7_self_contact/distance/PointLineDistance.py

@@ -1,3 +1,4 @@
+# ANCHOR: PL_val_grad
 import numpy as np
 
 def val(p, e0, e1):
@@ -23,6 +24,7 @@ def grad(p, e0, e1):
     g[4] = t27 + t24 * (p[1] - e0[1]) * 2.0
     g[5] = t26 - t24 * (p[0] - e0[0]) * 2.0
     return g
+# ANCHOR_END: PL_val_grad
 
 def hess(p, e0, e1):
     H = np.array([0.0] * 36)

7_self_contact/simulator.py

@@ -7,6 +7,7 @@
 import square_mesh   # square mesh
 import time_integrator
 
+# ANCHOR: sim_setup
 # simulation setup
 side_len = 0.45
 rho = 1000      # density of square
@@ -26,6 +27,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
+# ANCHOR_END: sim_setup
 [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

7_self_contact/square_mesh.py

@@ -23,6 +23,7 @@ def generate(side_length, n_seg):
 
     return [x, e]
 
+# ANCHOR: find_boundary
 def find_boundary(e):
     # index all half-edges for fast query
     edge_set = set()
@@ -41,6 +42,7 @@ def find_boundary(e):
             bp_set.add(eI[0])
             bp_set.add(eI[1])
     return [list(bp_set), be]
+# ANCHOR_END: find_boundary
 
 def write_to_file(frameNum, x, e):
     # Check if 'output' directory exists; if not, create it