ceil

Author: liminchen

mov_dirichlet/BarrierEnergy.py

@@ -6,24 +6,43 @@
 
 def val(x, n, o, contact_area):
     sum = 0.0
+    # floor:
     for i in range(0, len(x)):
         d = n.dot(x[i] - o)
         if d < dhat:
             s = d / dhat
             sum += contact_area[i] * dhat * kappa / 2 * (s - 1) * math.log(s)
+    # ceil:
+    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)
     return sum
 
 def grad(x, n, o, contact_area):
     g = np.array([[0.0, 0.0]] * len(x))
+    # floor:
     for i in range(0, len(x)):
         d = n.dot(x[i] - o)
         if d < dhat:
             s = d / dhat
             g[i] = contact_area[i] * dhat * (kappa / 2 * (math.log(s) / dhat + (s - 1) / d)) * n
+    # ceil:
+    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
+            local_grad = contact_area[i] * dhat * (kappa / 2 * (math.log(s) / dhat + (s - 1) / d)) * n
+            g[i] += local_grad
+            g[-1] -= local_grad
     return g
 
 def hess(x, n, o, contact_area):
     IJV = [[0] * 0, [0] * 0, np.array([0.0] * 0)]
+    # floor:
     for i in range(0, len(x)):
         d = n.dot(x[i] - o)
         if d < dhat:
@@ -33,14 +52,35 @@ 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:
+    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:
+            local_hess = contact_area[i] * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * np.outer(n, n)
+            index = [i, len(x) - 1]
+            for nI in range(0, 2):
+                for nJ in range(0, 2):
+                    for c in range(0, 2):
+                        for r in range(0, 2):
+                            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])
     return IJV
 
 def init_step_size(x, n, o, p):
     alpha = 1
+    # floor:
     for i in range(0, len(x)):
         p_n = p[i].dot(n)
         if p_n < 0:
             alpha = min(alpha, 0.9 * n.dot(x[i] - o) / -p_n)
+    # ceil:
+    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)
     return alpha
 
 def compute_mu_lambda(x, n, o, contact_area, mu):

mov_dirichlet/simulator.py

@@ -14,14 +14,15 @@
 n_seg = 4       # num of segments per side of the square
 h = 0.01        # time step size in s
 DBC = []        # no nodes need to be fixed
-ground_n = np.array([0.1, 1.0])     # normal of the slope
+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
-[x, e] = square_mesh.generate(side_len, n_seg)  # node positions and edge node indices
-v = np.array([[0.0, 0.0]] * len(x))             # velocity
+[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])
+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
 l2 = []
@@ -56,10 +57,13 @@ def screen_projection(x):
     # fill the background and draw the square
     screen.fill((255, 255, 255))
     pygame.draw.aaline(screen, (0, 0, 255), screen_projection([ground_o[0] - 3.0 * ground_n[1], ground_o[1] + 3.0 * ground_n[0]]), 
-        screen_projection([ground_o[0] + 3.0 * ground_n[1], ground_o[1] - 3.0 * ground_n[0]]))   # slope
+        screen_projection([ground_o[0] + 3.0 * ground_n[1], ground_o[1] - 3.0 * ground_n[0]]))   # ground
+    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:
         pygame.draw.aaline(screen, (0, 0, 255), screen_projection(x[eI[0]]), screen_projection(x[eI[1]]))
-    for xI in x:
+    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)
 
     pygame.display.flip()   # flip the display