added anchors

Author: liminchen

5_mov_dirichlet/BarrierEnergy.py

@@ -12,13 +12,14 @@ def val(x, n, o, contact_area):
         if d < dhat:
             s = d / dhat
             sum += contact_area[i] * dhat * kappa / 2 * (s - 1) * math.log(s)
-    # ceil:
+    # ANCHOR: ceiling_val
     n = np.array([0.0, -1.0])
     for i in range(0, len(x) - 1):
         d = n.dot(x[i] - x[-1])
         if d < dhat:
             s = d / dhat
             sum += contact_area[i] * dhat * kappa / 2 * (s - 1) * math.log(s)
+    # ANCHOR_END: ceiling_val
     return sum
 
 def grad(x, n, o, contact_area):
@@ -29,7 +30,7 @@ def grad(x, n, o, contact_area):
         if d < dhat:
             s = d / dhat
             g[i] = contact_area[i] * dhat * (kappa / 2 * (math.log(s) / dhat + (s - 1) / d)) * n
-    # ceil:
+    # ANCHOR: ceiling_grad
     n = np.array([0.0, -1.0])
     for i in range(0, len(x) - 1):
         d = n.dot(x[i] - x[-1])
@@ -38,6 +39,7 @@ def grad(x, n, o, 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_END: ceiling_grad
     return g
 
 def hess(x, n, o, contact_area):
@@ -52,7 +54,7 @@ def hess(x, n, o, contact_area):
                     IJV[0].append(i * 2 + r)
                     IJV[1].append(i * 2 + c)
                     IJV[2] = np.append(IJV[2], local_hess[r, c])
-    # ceil:
+    # ANCHOR: ceiling_hess
     n = np.array([0.0, -1.0])
     for i in range(0, len(x) - 1):
         d = n.dot(x[i] - x[-1])
@@ -66,6 +68,7 @@ def hess(x, n, o, 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_END: ceiling_hess
     return IJV
 
 def init_step_size(x, n, o, p):
@@ -75,12 +78,13 @@ def init_step_size(x, n, o, p):
         p_n = p[i].dot(n)
         if p_n < 0:
             alpha = min(alpha, 0.9 * n.dot(x[i] - o) / -p_n)
-    # ceil:
+    # ANCHOR: ceiling_ccd
     n = np.array([0.0, -1.0])
     for i in range(0, len(x) - 1):
         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_END: ceiling_ccd
     return alpha
 
 def compute_mu_lambda(x, n, o, contact_area, mu):

5_mov_dirichlet/simulator.py

@@ -13,17 +13,21 @@
 k = 2e4         # spring stiffness
 n_seg = 4       # num of segments per side of the square
 h = 0.01        # time step size in s
+# ANCHOR: ceiling_dbc_setup
 DBC = [(n_seg + 1) * (n_seg + 1)]       # dirichlet node index
 DBC_v = [np.array([0.0, -0.5])]         # dirichlet node velocity
 DBC_limit = [np.array([0.0, -0.6])]     # dirichlet node limit position
+# ANCHOR_END: ceiling_dbc_setup
 ground_n = np.array([0.0, 1.0])         # normal of the slope
 ground_n /= np.linalg.norm(ground_n)    # normalize ground normal vector just in case
 ground_o = np.array([0.0, -1.0])        # a point on the slope  
 mu = 0.11        # friction coefficient of the slope
 
 # initialize simulation
+# ANCHOR: ceiling_dof
 [x, e] = square_mesh.generate(side_len, n_seg)      # node positions and edge node indices
 x = np.append(x, [[0.0, side_len * 0.6]], axis=0)   # ceil origin (with normal [0.0, -1.0])
+# ANCHOR_END: ceiling_dof
 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
 # rest length squared

5_mov_dirichlet/time_integrator.py

@@ -17,17 +17,20 @@ def step_forward(x, e, v, m, l2, k, n, o, contact_area, mu, is_DBC, DBC, DBC_v,
     x_tilde = x + v * h     # implicit Euler predictive position
     x_n = copy.deepcopy(x)
     mu_lambda = BarrierEnergy.compute_mu_lambda(x, n, o, contact_area, mu)  # compute mu * lambda for each node using x^n
+    # ANCHOR: dbc_initialization
     DBC_target = [] # target position of each DBC in the current time step
     for i in range(0, len(DBC)):
         if (DBC_limit[i] - x_n[DBC[i]]).dot(DBC_v[i]) > 0:
             DBC_target.append(x_n[DBC[i]] + h * DBC_v[i])
         else:
             DBC_target.append(x_n[DBC[i]])
     DBC_stiff = 10  # initialize stiffness for DBC springs
+    # ANCHOR_END: dbc_initialization
 
     # Newton loop
     iter = 0
     E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
+    # ANCHOR: convergence_criteria
     [p, DBC_satisfied] = search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h)
     while (LA.norm(p, inf) / h > tol) | (sum(DBC_satisfied) != len(DBC)):   # also check whether all DBCs are satisfied
         print('Iteration', iter, ':')
@@ -37,6 +40,7 @@ def step_forward(x, e, v, m, l2, k, n, o, contact_area, mu, is_DBC, DBC, DBC_v,
             # increase DBC stiffness and recompute energy value record
             DBC_stiff *= 2
             E_last = IP_val(x, e, x_tilde, m, l2, k, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
+        # ANCHOR_END: convergence_criteria
 
         # filter line search
         alpha = BarrierEnergy.init_step_size(x, n, o, p)  # avoid interpenetration and tunneling
@@ -87,16 +91,20 @@ def IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_
 def search_dir(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, is_DBC, DBC, DBC_target, DBC_stiff, tol, h):
     projected_hess = IP_hess(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h)
     reshaped_grad = IP_grad(x, e, x_tilde, m, l2, k, n, o, contact_area, v, mu_lambda, DBC, DBC_target, DBC_stiff, h).reshape(len(x) * 2, 1)
+    # ANCHOR: dbc_check
     # check whether each DBC is satisfied
     DBC_satisfied = [False] * len(x)
     for i in range(0, len(DBC)):
         if LA.norm(x[DBC[i]] - DBC_target[i]) / h < tol:
             DBC_satisfied[DBC[i]] = True
+    # ANCHOR_END: dbc_check
+    # ANCHOR: dof_elimination
     # eliminate DOF if it's a satisfied DBC by modifying gradient and Hessian for DBC:
     for i, j in zip(*projected_hess.nonzero()):
         if (is_DBC[int(i / 2)] & DBC_satisfied[int(i / 2)]) | (is_DBC[int(j / 2)] & DBC_satisfied[int(i / 2)]): 
             projected_hess[i, j] = (i == j)
     for i in range(0, len(x)):
         if is_DBC[i] & DBC_satisfied[i]:
             reshaped_grad[i * 2] = reshaped_grad[i * 2 + 1] = 0.0
-    return [spsolve(projected_hess, -reshaped_grad).reshape(len(x), 2), DBC_satisfied]
+    return [spsolve(projected_hess, -reshaped_grad).reshape(len(x), 2), DBC_satisfied]
+    # ANCHOR_END: dof_elimination