Vectorized 'do_initial_projection()'.

scripts/tcd/DG_Approximation.py

@@ -181,26 +181,25 @@ def do_initial_projection(init_cond, mesh, basis, quadrature,
         number of mesh cells and p being the polynomial degree of the basis.
 
     """
-    # Initialize matrix and set first entry to accommodate for ghost cell
-    output_matrix = [0]
-
-    for eval_point in mesh.non_ghost_cells:
-        new_row = []
-        for degree in range(basis.polynomial_degree + 1):
-            new_row.append(np.float64(sum(init_cond.calculate(
-                x=eval_point + mesh.cell_len/2
-                * quadrature.nodes[point] - x_shift, mesh=mesh)
-                * basis.basis[degree].subs(
-                    x, quadrature.nodes[point])
-                * quadrature.weights[point]
-                for point in range(quadrature.num_nodes))))
-
-        new_row = np.array(new_row)
-        output_matrix.append(basis.inverse_mass_matrix @ new_row)
-
-    # Set ghost cells to respective value
-    output_matrix[0] = output_matrix[mesh.num_cells]
-    output_matrix.append(output_matrix[1])
-
-    # print(np.array(output_matrix).shape)
-    return np.transpose(np.array(output_matrix))
+    # Initialize output matrix
+    output_matrix = np.zeros((basis.polynomial_degree+1, mesh.num_cells+2))
+
+    # Calculate basis based on quadrature
+    basis_matrix = np.array([np.vectorize(basis.basis[degree].subs)(
+        x, quadrature.nodes) for degree in range(basis.polynomial_degree+1)])
+
+    # Calculate points based on initial condition
+    points = quadrature.nodes[:, np.newaxis] @ \
+        np.repeat(mesh.cell_len/2, mesh.num_cells)[:, np.newaxis].T + \
+        np.tile(mesh.non_ghost_cells - x_shift, (quadrature.num_nodes, 1))
+    init_matrix = np.vectorize(init_cond.calculate, otypes=[np.float])(
+        x=points, mesh=mesh)
+
+    # Set output matrix for regular cells
+    output_matrix[:, 1:-1] = (basis_matrix * quadrature.weights) @ init_matrix
+
+    # Set ghost cells
+    output_matrix[:, 0] = output_matrix[:, -2]
+    output_matrix[:, -1] = output_matrix[:, 1]
+
+    return output_matrix