subspace simulation

Author: liminchen

9_reduced_DOF/simulator.py

@@ -4,6 +4,7 @@
 import pygame       # pygame for visualization
 pygame.init()
 
+import utils
 import square_mesh   # square mesh
 import time_integrator
 
@@ -33,6 +34,8 @@
 # ANCHOR: contact_area
 contact_area = [side_len / n_seg] * len(x)     # perimeter split to each node
 # ANCHOR_END: contact_area
+# compute reduced basis using 0: polynomial functions or 1: modal reduction
+reduced_basis = utils.compute_reduced_basis(x, e, l2, k, is_DBC, method=0, order=3)
 
 # simulation with visualization
 resolution = np.array([900, 900])
@@ -64,7 +67,7 @@ def screen_projection(x):
     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, y_ground, contact_area, is_DBC, h, 1e-2)
+    [x, v] = time_integrator.step_forward(x, e, v, m, l2, k, y_ground, contact_area, is_DBC, reduced_basis, h, 1e-2)
     time_step += 1
     pygame.time.wait(int(h * 1000))
     square_mesh.write_to_file(time_step, x, n_seg)

9_reduced_DOF/time_integrator.py

@@ -11,14 +11,14 @@
 import GravityEnergy
 import BarrierEnergy
 
-def step_forward(x, e, v, m, l2, k, y_ground, contact_area, is_DBC, h, tol):
+def step_forward(x, e, v, m, l2, k, y_ground, contact_area, is_DBC, reduced_basis, h, tol):
     x_tilde = x + v * h     # implicit Euler predictive position
     x_n = copy.deepcopy(x)
 
     # Newton loop
     iter = 0
     E_last = IP_val(x, e, x_tilde, m, l2, k, y_ground, contact_area, h)
-    p = search_dir(x, e, x_tilde, m, l2, k, y_ground, contact_area, is_DBC, h)
+    p = search_dir(x, e, x_tilde, m, l2, k, y_ground, contact_area, is_DBC, reduced_basis, h)
     while LA.norm(p, inf) / h > tol:
         print('Iteration', iter, ':')
         print('residual =', LA.norm(p, inf) / h)
@@ -33,7 +33,7 @@ def step_forward(x, e, v, m, l2, k, y_ground, contact_area, is_DBC, h, tol):
 
         x += alpha * p
         E_last = IP_val(x, e, x_tilde, m, l2, k, y_ground, contact_area, h)
-        p = search_dir(x, e, x_tilde, m, l2, k, y_ground, contact_area, is_DBC, h)
+        p = search_dir(x, e, x_tilde, m, l2, k, y_ground, contact_area, is_DBC, reduced_basis, h)
         iter += 1
 
     v = (x - x_n) / h   # implicit Euler velocity update
@@ -56,7 +56,7 @@ def IP_hess(x, e, x_tilde, m, l2, k, y_ground, contact_area, h):
     H = sparse.coo_matrix((IJV[2], (IJV[0], IJV[1])), shape=(len(x) * 2, len(x) * 2)).tocsr()
     return H
 
-def search_dir(x, e, x_tilde, m, l2, k, y_ground, contact_area, is_DBC, h):
+def search_dir(x, e, x_tilde, m, l2, k, y_ground, contact_area, is_DBC, reduced_basis, h):
     projected_hess = IP_hess(x, e, x_tilde, m, l2, k, y_ground, contact_area, h)
     reshaped_grad = IP_grad(x, e, x_tilde, m, l2, k, y_ground, contact_area, h).reshape(len(x) * 2, 1)
     # eliminate DOF by modifying gradient and Hessian for DBC:
@@ -66,4 +66,6 @@ def search_dir(x, e, x_tilde, m, l2, k, y_ground, contact_area, is_DBC, h):
     for i in range(0, len(x)):
         if is_DBC[i]:
             reshaped_grad[i * 2] = reshaped_grad[i * 2 + 1] = 0.0
-    return spsolve(projected_hess, -reshaped_grad).reshape(len(x), 2)
+    reduced_hess = reduced_basis.T.dot(projected_hess.dot(reduced_basis)) # applying chain rule
+    reduced_grad = reduced_basis.T.dot(reshaped_grad) # applying chain rule
+    return (reduced_basis.dot(spsolve(reduced_hess, -reduced_grad))).reshape(len(x), 2) # transform to full space after the solve

9_reduced_DOF/utils.py

@@ -1,9 +1,63 @@
 import numpy as np
 import numpy.linalg as LA
 
+import scipy.sparse as sparse
+from scipy.sparse.linalg import eigsh
+
+import MassSpringEnergy
+
 def make_PSD(hess):
     [lam, V] = LA.eigh(hess)    # Eigen decomposition on symmetric matrix
     # set all negative Eigenvalues to 0
     for i in range(0, len(lam)):
         lam[i] = max(0, lam[i])
-    return np.matmul(np.matmul(V, np.diag(lam)), np.transpose(V))
+    return np.matmul(np.matmul(V, np.diag(lam)), np.transpose(V))
+
+def compute_reduced_basis(x, e, l2, k, is_DBC, method, order):
+    if method == 0: # polynomial basis
+        if order == 1: # linear basis, or affine basis
+            basis = np.zeros((len(x) * 2, 6)) # 1, x, y for both x- and y-displacements
+            for i in range(len(x)):
+                for d in range(2):
+                    if not is_DBC[i]: # ignore the floor DOF
+                        basis[i * 2 + d][d * 3] = 1
+                        basis[i * 2 + d][d * 3 + 1] = x[i][0]
+                        basis[i * 2 + d][d * 3 + 2] = x[i][1]
+        elif order == 2: # quadratic polynomial basis 
+            basis = np.zeros((len(x) * 2, 12)) # 1, x, y, x^2, xy, y^2 for both x- and y-displacements
+            for i in range(len(x)):
+                for d in range(2):
+                    if not is_DBC[i]: # ignore the floor DOF
+                        basis[i * 2 + d][d * 6] = 1
+                        basis[i * 2 + d][d * 6 + 1] = x[i][0]
+                        basis[i * 2 + d][d * 6 + 2] = x[i][1]
+                        basis[i * 2 + d][d * 6 + 3] = x[i][0] * x[i][0]
+                        basis[i * 2 + d][d * 6 + 4] = x[i][0] * x[i][1]
+                        basis[i * 2 + d][d * 6 + 5] = x[i][1] * x[i][1]
+        elif order == 3: # cubic polynomial basis
+            basis = np.zeros((len(x) * 2, 20)) # 1, x, y, x^2, xy, y^2, x^3, x^2y, xy^2, y^3 for both x- and y-displacements
+            for i in range(len(x)):
+                for d in range(2):
+                    if not is_DBC[i]: # ignore the floor DOF
+                        basis[i * 2 + d][d * 10] = 1
+                        basis[i * 2 + d][d * 10 + 1] = x[i][0]
+                        basis[i * 2 + d][d * 10 + 2] = x[i][1]
+                        basis[i * 2 + d][d * 10 + 3] = x[i][0] * x[i][0]
+                        basis[i * 2 + d][d * 10 + 4] = x[i][0] * x[i][1]
+                        basis[i * 2 + d][d * 10 + 5] = x[i][1] * x[i][1]
+                        basis[i * 2 + d][d * 10 + 6] = x[i][0] * x[i][0] * x[i][0]
+                        basis[i * 2 + d][d * 10 + 7] = x[i][0] * x[i][0] * x[i][1]
+                        basis[i * 2 + d][d * 10 + 8] = x[i][0] * x[i][1] * x[i][1]
+                        basis[i * 2 + d][d * 10 + 9] = x[i][1] * x[i][1] * x[i][1]
+        else:
+            print("unsupported order of polynomial basis for reduced DOF")
+            exit()
+        return basis
+    else: # modal-order reduction
+        IJV = MassSpringEnergy.hess(x, e, l2, k) #TODO: no SPD projection, switch to NH
+        H = sparse.coo_matrix((IJV[2], (IJV[0], IJV[1])), shape=(len(x) * 2, len(x) * 2)).tocsr()
+        for i, j in zip(*H.nonzero()):
+            if is_DBC[int(i / 2)] | is_DBC[int(j / 2)]: 
+                H[i, j] = (i == j) # ignore the floor DOF
+        eigenvalues, eigenvectors = eigsh(H, k=order, which='SM') # get 'order' eigenvectors with smallest eigenvalues 
+        return eigenvectors