From abc8e917d57b416ec49b4904eb1441bd0143cfef Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?K=C3=BChle=2C=20Laura=20Christine=20=28lakue103=29?=
 <laura.kuehle@uni-duesseldorf.de>
Date: Wed, 9 Nov 2022 18:38:12 +0100
Subject: [PATCH] Vectorized '_reset()' for update scheme.

---
 Snakefile                    |  2 +-
 scripts/tcd/Update_Scheme.py | 45 +++++++++++++++---------------------
 2 files changed, 19 insertions(+), 28 deletions(-)

diff --git a/Snakefile b/Snakefile
index 10e4ed1..7bcb15c 100644
--- a/Snakefile
+++ b/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
diff --git a/scripts/tcd/Update_Scheme.py b/scripts/tcd/Update_Scheme.py
index 771ecb7..b418732 100644
--- a/scripts/tcd/Update_Scheme.py
+++ b/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):
-                if self._wave_speed > 0:
-                    new_entry = -1.0
-                    if (j < i) & ((i+j) % 2 == 1):
-                        new_entry = 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.ones(matrix_shape)
+        if self._wave_speed > 0:
+            matrix[np.fromfunction(
+                lambda i, j: (j >= i) | ((i+j) % 2 == 0), matrix_shape)] = -1.0
+        else:
+            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 \
-                    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)
+        matrix = np.fromfunction(lambda i, j: (-1.0)**i, matrix_shape) \
+            if self._wave_speed > 0 \
+            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."""
-- 
GitLab