Vectorized '_reset()' for update scheme.

Snakefile

@@ -23,7 +23,7 @@ TODO: Discuss name for quadrature mesh (now: grid)
 TODO: Contemplate using lambdify for basis
 
 Urgent:
-TODO: Vectorize '_reset()' for update scheme
+TODO: Vectorize '_reset()' for update scheme -> Done
 TODO: Vectorize 'plot_details()'
 TODO: Replace loops/list comprehension with vectorization if feasible
 TODO: Replace loops with list comprehension if feasible

scripts/tcd/Update_Scheme.py

@@ -55,36 +55,27 @@ class UpdateScheme(ABC):
 
     def _reset(self):
         """Resets instance variables."""
+        matrix_shape = (self._polynomial_degree+1, self._polynomial_degree+1)
+        root_vector = np.array([np.sqrt(degree+0.5) for degree in range(
+            self._polynomial_degree+1)])
+        degree_matrix = np.matmul(root_vector[:, np.newaxis],
+                                  root_vector[:, np.newaxis].T)
+
         # Set stiffness matrix
-        matrix = []
-        for i in range(self._polynomial_degree+1):
-            new_row = []
-            for j in range(self._polynomial_degree+1):
+        matrix = np.ones(matrix_shape)
         if self._wave_speed > 0:
-                    new_entry = -1.0
-                    if (j < i) & ((i+j) % 2 == 1):
-                        new_entry = 1.0
+            matrix[np.fromfunction(
+                lambda i, j: (j >= i) | ((i+j) % 2 == 0), matrix_shape)] = -1.0
         else:
-                    new_entry = 1.0
-                    if (j > i) & ((i+j) % 2 == 1):
-                        new_entry = -1.0
-                new_row.append(new_entry*np.sqrt((i+0.5) * (j+0.5)))
-            matrix.append(new_row)
-        self._stiffness_matrix = np.array(matrix)
-        # former: inv_mass @ np.array(matrix)
+            matrix[np.fromfunction(
+                lambda i, j: (j > i) & ((i+j) % 2 == 1), matrix_shape)] = -1.0
+        self._stiffness_matrix = matrix * degree_matrix
 
         # Set boundary matrix
-        matrix = []
-        for i in range(self._polynomial_degree+1):
-            new_row = []
-            for j in range(self._polynomial_degree+1):
-                new_entry = np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**i \
+        matrix = np.fromfunction(lambda i, j: (-1.0)**i, matrix_shape) \
             if self._wave_speed > 0 \
-                    else np.sqrt((i+0.5) * (j+0.5)) * (-1.0)**j
-                new_row.append(new_entry)
-            matrix.append(new_row)
-        self._boundary_matrix = np.array(matrix)
-        # former: inv_mass @ np.array(matrix)
+            else np.fromfunction(lambda i, j: (-1.0)**j, matrix_shape)
+        self._boundary_matrix = matrix * degree_matrix
 
     def get_name(self):
         """Returns string of class name."""