added anchors

Author: liminchen

6_inv_free/NeoHookeanEnergy.py

@@ -1,3 +1,4 @@
+# ANCHOR: helper_func
 import utils
 import numpy as np
 import math
@@ -78,7 +79,9 @@ def d2Psi_div_dF2(F, mu, lam):
                         + M[3, 0] * U[i, 1] * VT[1, j] * U[r, 0] * VT[0, s] \
                         + M[3, 3] * U[i, 1] * VT[1, j] * U[r, 1] * VT[1, s]
     return dP_div_dF
+# ANCHOR_END: helper_func
 
+# ANCHOR: stress_deriv
 def deformation_grad(x, elemVInd, IB):
     F = [x[elemVInd[1]] - x[elemVInd[0]], x[elemVInd[2]] - x[elemVInd[0]]]
     return np.transpose(F).dot(IB)
@@ -124,7 +127,9 @@ def d2Psi_div_dx2(dP_div_dF, IB):  # applying chain-rule, d2Psi_div_dx2 = dF_div
         result[0, colI] = -_000 - _101 - _010 - _111
         result[1, colI] = -_200 - _301 - _210 - _311
     return result
+# ANCHOR_END: stress_deriv
 
+# ANCHOR: val_grad_hess
 def val(x, e, vol, IB, mu, lam):
     sum = 0.0
     for i in range(0, len(e)):
@@ -157,7 +162,9 @@ def hess(x, e, vol, IB, mu, lam):
                         IJV[1][ind] = e[i][xJ] * 2 + dJ
                         IJV[2][ind] = local_hess[xI * 2 + dI, xJ * 2 + dJ]
     return IJV
+# ANCHOR_END: val_grad_hess
 
+# ANCHOR: filter_line_search
 def init_step_size(x, e, p):
     alpha = 1
     for i in range(0, len(e)):
@@ -172,4 +179,5 @@ def init_step_size(x, e, p):
         critical_alpha = utils.smallest_positive_real_root_quad(a, b, c)
         if critical_alpha > 0:
             alpha = min(alpha, critical_alpha)
-    return alpha
+    return alpha
+# ANCHOR_END: filter_line_search

6_inv_free/simulator.py

@@ -10,8 +10,10 @@
 # simulation setup
 side_len = 1
 rho = 1000      # density of square
+# ANCHOR: lame_param
 E = 1e5         # Young's modulus
 nu = 0.4        # Poisson's ratio
+# ANCHOR_END: lame_param
 n_seg = 4       # num of segments per side of the square
 h = 0.01        # time step size in s
 DBC = [(n_seg + 1) * (n_seg + 1)]       # dirichlet node index
@@ -27,6 +29,7 @@
 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
+# ANCHOR: elem_precomp
 # rest shape basis, volume, and lame parameters
 vol = [0.0] * len(e)
 IB = [np.array([[0.0, 0.0]] * 2)] * len(e)
@@ -36,6 +39,7 @@
     IB[i] = np.linalg.inv(np.transpose(TB))
 mu_lame = [0.5 * E / (1 + nu)] * len(e)
 lam = [E * nu / ((1 + nu) * (1 - 2 * nu))] * len(e)
+# ANCHOR_END: elem_precomp
 # identify whether a node is Dirichlet
 is_DBC = [False] * len(x)
 for i in DBC:
@@ -68,9 +72,11 @@ def screen_projection(x):
     pygame.draw.aaline(screen, (0, 0, 255), screen_projection([x[-1][0] + 3.0, x[-1][1]]), 
         screen_projection([x[-1][0] - 3.0, x[-1][1]]))   # ceil
     for eI in e:
+        # ANCHOR: draw_tri
         pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[0]]), screen_projection(x[eI[1]]))
         pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[1]]), screen_projection(x[eI[2]]))
         pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[2]]), screen_projection(x[eI[0]]))
+        # ANCHOR_END: draw_tri
     for xId in range(0, len(x) - 1):
         xI = x[xId]
         pygame.draw.circle(screen, (0, 0, 255), screen_projection(xI), 0.1 * side_len / n_seg * scale)

6_inv_free/square_mesh.py

@@ -9,6 +9,7 @@ def generate(side_length, n_seg):
         for j in range(0, n_seg + 1):
             x[i * (n_seg + 1) + j] = [-side_length / 2 + i * step, -side_length / 2 + j * step]
 
+    # ANCHOR: tri_vert_ind
     # connect the nodes with triangle elements
     e = []
     for i in range(0, n_seg):
@@ -20,6 +21,7 @@ def generate(side_length, n_seg):
             else:
                 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])
+    # ANCHOR_END: tri_vert_ind
 
     return [x, e]
 

6_inv_free/time_integrator.py

@@ -39,7 +39,9 @@ def step_forward(x, e, v, m, vol, IB, mu_lame, lam, n, o, contact_area, mu, is_D
             E_last = IP_val(x, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h)
 
         # filter line search
+        # ANCHOR: apply_filter
         alpha = min(BarrierEnergy.init_step_size(x, n, o, p), NeoHookeanEnergy.init_step_size(x, e, p))  # avoid interpenetration, tunneling, and inversion
+        # ANCHOR_END: apply_filter
         while IP_val(x + alpha * p, e, x_tilde, m, vol, IB, mu_lame, lam, n, o, contact_area, (x + alpha * p - x_n) / h, mu_lambda, DBC, DBC_target, DBC_stiff, h) > E_last:
             alpha /= 2
         print('step size =', alpha)

6_inv_free/utils.py

@@ -9,6 +9,7 @@ def make_PSD(hess):
         lam[i] = max(0, lam[i])
     return np.matmul(np.matmul(V, np.diag(lam)), np.transpose(V))
 
+# ANCHOR: find_positive_real_root
 def smallest_positive_real_root_quad(a, b, c, tol = 1e-6):
     # return negative value if no positive real root is found
     t = 0
@@ -25,4 +26,5 @@ def smallest_positive_real_root_quad(a, b, c, tol = 1e-6):
                 t = (-b + math.sqrt(desc)) / (2 * a)
         else: # desv<0 ==> imag, f(x) > 0 for all x > 0
             t = -1
-    return t
+    return t
+# ANCHOR_END: find_positive_real_root