init inv free

Author: liminchen

inv_free/BarrierEnergy.py

@@ -0,0 +1,93 @@
+import math
+import numpy as np
+
+dhat = 0.01
+kappa = 1e5
+
+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:
+            local_hess = contact_area[i] * dhat * kappa / (2 * d * d * dhat) * (d + dhat) * np.outer(n, n)
+            for c in range(0, 2):
+                for r in range(0, 2):
+                    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):
+    mu_lambda = np.array([0.0] * len(x))
+    for i in range(0, len(x)):
+        d = n.dot(x[i] - o)
+        if d < dhat:
+            s = d / dhat
+            mu_lambda[i] = mu * -contact_area[i] * dhat * (kappa / 2 * (math.log(s) / dhat + (s - 1) / d))
+    return mu_lambda

inv_free/FrictionEnergy.py

@@ -0,0 +1,60 @@
+import numpy as np
+import utils
+
+epsv = 1e-3
+
+def f0(vbarnorm, epsv, hhat):
+    if vbarnorm >= epsv:
+        return vbarnorm * hhat
+    else:
+        vbarnormhhat = vbarnorm * hhat
+        epsvhhat = epsv * hhat
+        return vbarnormhhat * vbarnormhhat * (-vbarnormhhat / 3.0 + epsvhhat) / (epsvhhat * epsvhhat) + epsvhhat / 3.0
+
+def f1_div_vbarnorm(vbarnorm, epsv):
+    if vbarnorm >= epsv:
+        return 1.0 / vbarnorm
+    else:
+        return (-vbarnorm + 2.0 * epsv) / (epsv * epsv)
+
+def f_hess_term(vbarnorm, epsv):
+    if vbarnorm >= epsv:
+        return -1.0 / (vbarnorm * vbarnorm)
+    else:
+        return -1.0 / (epsv * epsv)
+
+def val(v, mu_lambda, hhat, n):
+    sum = 0.0
+    T = np.identity(2) - np.outer(n, n) # tangent of slope is constant
+    for i in range(0, len(v)):
+        if mu_lambda[i] > 0:
+            vbar = np.transpose(T).dot(v[i])
+            sum += mu_lambda[i] * f0(np.linalg.norm(vbar), epsv, hhat)
+    return sum
+
+def grad(v, mu_lambda, hhat, n):
+    g = np.array([[0.0, 0.0]] * len(v))
+    T = np.identity(2) - np.outer(n, n) # tangent of slope is constant
+    for i in range(0, len(v)):
+        if mu_lambda[i] > 0:
+            vbar = np.transpose(T).dot(v[i])
+            g[i] = mu_lambda[i] * f1_div_vbarnorm(np.linalg.norm(vbar), epsv) * T.dot(vbar)
+    return g
+
+def hess(v, mu_lambda, hhat, n):
+    IJV = [[0] * 0, [0] * 0, np.array([0.0] * 0)]
+    T = np.identity(2) - np.outer(n, n) # tangent of slope is constant
+    for i in range(0, len(v)):
+        if mu_lambda[i] > 0:
+            vbar = np.transpose(T).dot(v[i])
+            vbarnorm = np.linalg.norm(vbar)
+            inner_term = f1_div_vbarnorm(vbarnorm, epsv) * np.identity(2)
+            if vbarnorm != 0:
+                inner_term += f_hess_term(vbarnorm, epsv) / vbarnorm * np.outer(vbar, vbar)
+            local_hess = mu_lambda[i] * T.dot(utils.make_PD(inner_term)).dot(np.transpose(T)) / hhat
+            for c in range(0, 2):
+                for r in range(0, 2):
+                    IJV[0].append(i * 2 + r)
+                    IJV[1].append(i * 2 + c)
+                    IJV[2] = np.append(IJV[2], local_hess[r, c])
+    return IJV

inv_free/GravityEnergy.py

@@ -0,0 +1,17 @@
+import numpy as np
+
+gravity = [0.0, -9.81]
+
+def val(x, m):
+    sum = 0.0
+    for i in range(0, len(x)):
+        sum += -m[i] * x[i].dot(gravity)
+    return sum
+
+def grad(x, m):
+    g = np.array([gravity] * len(x))
+    for i in range(0, len(x)):
+        g[i] *= -m[i]
+    return g
+
+# Hessian is 0

inv_free/InertiaEnergy.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+def val(x, x_tilde, m):
+    sum = 0.0
+    for i in range(0, len(x)):
+        diff = x[i] - x_tilde[i]
+        sum += 0.5 * m[i] * diff.dot(diff)
+    return sum
+
+def grad(x, x_tilde, m):
+    g = np.array([[0.0, 0.0]] * len(x))
+    for i in range(0, len(x)):
+        g[i] = m[i] * (x[i] - x_tilde[i])
+    return g
+
+def hess(x, x_tilde, m):
+    IJV = [[0] * (len(x) * 2), [0] * (len(x) * 2), np.array([0.0] * (len(x) * 2))]
+    for i in range(0, len(x)):
+        for d in range(0, 2):
+            IJV[0][i * 2 + d] = i * 2 + d
+            IJV[1][i * 2 + d] = i * 2 + d
+            IJV[2][i * 2 + d] = m[i]
+    return IJV

inv_free/MassSpringEnergy.py

@@ -0,0 +1,35 @@
+import numpy as np
+import utils
+
+def val(x, e, l2, k):
+    sum = 0.0
+    for i in range(0, len(e)):
+        diff = x[e[i][0]] - x[e[i][1]]
+        sum += l2[i] * 0.5 * k[i] * (diff.dot(diff) / l2[i] - 1) ** 2
+    return sum
+
+def grad(x, e, l2, k):
+    g = np.array([[0.0, 0.0]] * len(x))
+    for i in range(0, len(e)):
+        diff = x[e[i][0]] - x[e[i][1]]
+        g_diff = 2 * k[i] * (diff.dot(diff) / l2[i] - 1) * diff
+        g[e[i][0]] += g_diff
+        g[e[i][1]] -= g_diff
+    return g
+
+def hess(x, e, l2, k):
+    IJV = [[0] * (len(e) * 16), [0] * (len(e) * 16), np.array([0.0] * (len(e) * 16))]
+    for i in range(0, len(e)):
+        diff = x[e[i][0]] - x[e[i][1]]
+        H_diff = 2 * k[i] / l2[i] * (2 * np.outer(diff, diff) + (diff.dot(diff) - l2[i]) * np.identity(2))
+        H_local = utils.make_PD(np.block([[H_diff, -H_diff], [-H_diff, H_diff]]))
+        # add to global matrix
+        for nI in range(0, 2):
+            for nJ in range(0, 2):
+                indStart = i * 16 + (nI * 2 + nJ) * 4
+                for r in range(0, 2):
+                    for c in range(0, 2):
+                        IJV[0][indStart + r * 2 + c] = e[i][nI] * 2 + r
+                        IJV[1][indStart + r * 2 + c] = e[i][nJ] * 2 + c
+                        IJV[2][indStart + r * 2 + c] = H_local[nI * 2 + r, nJ * 2 + c]
+    return IJV

inv_free/SpringEnergy.py

@@ -0,0 +1,23 @@
+import numpy as np
+
+def val(x, m, DBC, DBC_target, k):
+    sum = 0.0
+    for i in range(0, len(DBC)):
+        diff = x[DBC[i]] - DBC_target[i]
+        sum += 0.5 * k * m[DBC[i]] * diff.dot(diff)
+    return sum
+
+def grad(x, m, DBC, DBC_target, k):
+    g = np.array([[0.0, 0.0]] * len(x))
+    for i in range(0, len(DBC)):
+        g[DBC[i]] = k * m[DBC[i]] * (x[DBC[i]] - DBC_target[i])
+    return g
+
+def hess(x, m, DBC, DBC_target, k):
+    IJV = [[0] * 0, [0] * 0, np.array([0.0] * 0)]
+    for i in range(0, len(DBC)):
+        for d in range(0, 2):
+            IJV[0].append(DBC[i] * 2 + d)
+            IJV[1].append(DBC[i] * 2 + d)
+            IJV[2] = np.append(IJV[2], k * m[DBC[i]])
+    return IJV

inv_free/readme.md

@@ -0,0 +1,14 @@
+# Inversion-free Hyperelastic Solids Simulation
+
+A square falling onto a ground under gravity and then compressed by a ceiling is simulated with an inversion-free hyperelastic potential and implicit Euler time integration.
+Each time step is solved by minimizing the Incremental Potential with the projected Newton method.
+
+## Dependencies
+```
+pip install numpy scipy pygame
+```
+
+## Run
+```
+python simulator.py
+```

inv_free/simulator.py

@@ -0,0 +1,78 @@
+# Mass-Spring Solids Simulation
+
+import numpy as np  # numpy for linear algebra
+import pygame       # pygame for visualization
+pygame.init()
+
+import square_mesh   # square mesh
+import time_integrator
+
+# simulation setup
+side_len = 1
+rho = 1000      # density of square
+k = 2e4         # spring stiffness
+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
+DBC_v = [np.array([0.0, -0.5])]         # dirichlet node velocity
+DBC_limit = [np.array([0.0, -0.6])]     # dirichlet node limit position
+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
+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 = []
+for i in range(0, len(e)):
+    diff = x[e[i][0]] - x[e[i][1]]
+    l2.append(diff.dot(diff))
+k = [k] * len(e)    # spring stiffness
+# identify whether a node is Dirichlet
+is_DBC = [False] * len(x)
+for i in DBC:
+    is_DBC[i] = True
+contact_area = [side_len / n_seg] * len(x)     # perimeter split to each node
+
+# simulation with visualization
+resolution = np.array([900, 900])
+offset = resolution / 2
+scale = 200
+def screen_projection(x):
+    return [offset[0] + scale * x[0], resolution[1] - (offset[1] + scale * x[1])]
+
+time_step = 0
+screen = pygame.display.set_mode(resolution)
+running = True
+while running:
+    # run until the user asks to quit
+    for event in pygame.event.get():
+        if event.type == pygame.QUIT:
+            running = False
+    
+    print('### Time step', time_step, '###')
+
+    # 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]]))   # 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 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
+
+    # step forward simulation and wait for screen refresh
+    [x, v] = time_integrator.step_forward(x, e, v, m, l2, k, ground_n, ground_o, contact_area, mu, is_DBC, DBC, DBC_v, DBC_limit, h, 1e-2)
+    time_step += 1
+    pygame.time.wait(int(h * 1000))
+
+pygame.quit()

inv_free/square_mesh.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+def generate(side_length, n_seg):
+    # sample nodes uniformly on a square
+    x = np.array([[0.0, 0.0]] * ((n_seg + 1) ** 2))
+    step = side_length / n_seg
+    for i in range(0, n_seg + 1):
+        for j in range(0, n_seg + 1):
+            x[i * (n_seg + 1) + j] = [-side_length / 2 + i * step, -side_length / 2 + j * step]
+    
+    # connect the nodes with edges
+    e = []
+    # horizontal edges
+    for i in range(0, n_seg):
+        for j in range(0, n_seg + 1):
+            e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j])
+    # vertical edges
+    for i in range(0, n_seg + 1):
+        for j in range(0, n_seg):
+            e.append([i * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
+    # diagonals
+    for i in range(0, n_seg):
+        for j in range(0, n_seg):
+            e.append([i * (n_seg + 1) + j, (i + 1) * (n_seg + 1) + j + 1])
+            e.append([(i + 1) * (n_seg + 1) + j, i * (n_seg + 1) + j + 1])
+
+    return [x, e]