diff --git a/scripts/tcd/DG_Approximation.py b/scripts/tcd/DG_Approximation.py index 0b64d9d6031ab66ecd6a707d0cbac6b025d40f87..a816f08df9df835f1752d498946b9440ac6eefd1 100644 --- a/scripts/tcd/DG_Approximation.py +++ b/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