diff --git a/scripts/tcd/ANN_Data_Generator.py b/scripts/tcd/ANN_Data_Generator.py
index 881f4d29e4fde4236af690c0277f8e07521d1f82..783f3665638778c9bd3ffd99e9fa572248437f4e 100644
--- a/scripts/tcd/ANN_Data_Generator.py
+++ b/scripts/tcd/ANN_Data_Generator.py
@@ -7,8 +7,8 @@ import os
 import time
 import numpy as np
 
-from .DG_Approximation import do_initial_projection
 from .Mesh import Mesh
+from .projection_utils import do_initial_projection
 from .Quadrature import Gauss
 from .Basis_Function import OrthonormalLegendre
 
diff --git a/scripts/tcd/DG_Approximation.py b/scripts/tcd/DG_Approximation.py
index 7df9b44e11f786350047513f4f3f221ebb4a9c8c..a4cf1df3f853e6af8aa3ca26f05bb7d8eb994ff2 100644
--- a/scripts/tcd/DG_Approximation.py
+++ b/scripts/tcd/DG_Approximation.py
@@ -4,20 +4,17 @@
 
 """
 import json
-import numpy as np
-from sympy import Symbol
 import math
 
 from .Basis_Function import Basis
 from .encoding_utils import encode_ndarray
 from .Initial_Condition import InitialCondition
 from .Mesh import Mesh
+from .projection_utils import do_initial_projection
 from .Quadrature import Quadrature
 from .Troubled_Cell_Detector import TroubledCellDetector
 from .Update_Scheme import UpdateScheme
 
-x = Symbol('x')
-
 
 class DGScheme:
     """Class for Discontinuous Galerkin Method.
@@ -152,59 +149,3 @@ class DGScheme:
         with open(data_file + '.json', 'w') \
                 as json_file:
             json_file.write(json.dumps(approx_stats))
-
-
-def do_initial_projection(init_cond, mesh, basis, quadrature,
-                          x_shift=0):
-    """Calculates initial projection.
-
-    Calculates a projection at time step 0 and adds ghost cells on both
-    sides of the array.
-
-    Parameters
-    ----------
-    init_cond : InitialCondition object
-        Initial condition used for calculation.
-    mesh : Mesh object
-        Mesh for calculation.
-    basis: Basis object
-        Basis used for calculation.
-    quadrature: Quadrature object
-        Quadrature used for evaluation.
-    x_shift: float, optional
-        Shift in x-direction for each evaluation point. Default: 0.
-
-    Returns
-    -------
-    ndarray
-        Matrix containing projection of size (N+2, p+1) with N being the
-        number of mesh cells and p being the polynomial degree of the basis.
-
-    """
-    # Initialize output matrix
-    output_matrix = np.zeros((basis.polynomial_degree+1,
-                              mesh.num_cells+2*mesh.num_ghost_cells))
-
-    # 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[:, mesh.num_ghost_cells:
-                  output_matrix.shape[-1]-mesh.num_ghost_cells] = \
-        (basis_matrix * quadrature.weights) @ init_matrix
-
-    # Set ghost cells
-    output_matrix[:, :mesh.num_ghost_cells] = \
-        output_matrix[:, -2*mesh.num_ghost_cells:-mesh.num_ghost_cells]
-    output_matrix[:, output_matrix.shape[-1]-mesh.num_ghost_cells:] = \
-        output_matrix[:, mesh.num_ghost_cells:2*mesh.num_ghost_cells]
-
-    return output_matrix
diff --git a/scripts/tcd/projection_utils.py b/scripts/tcd/projection_utils.py
index a7c061ea21d4e16efc2874b35d0e42f6aa6e3a08..17e0b3e1dcde2d33ac1bf0d86868fabc11d361d0 100644
--- a/scripts/tcd/projection_utils.py
+++ b/scripts/tcd/projection_utils.py
@@ -91,3 +91,59 @@ def calculate_exact_solution(
     exact = np.reshape(exact, (1, exact.size))
 
     return grid, exact
+
+
+def do_initial_projection(init_cond, mesh, basis, quadrature,
+                          x_shift=0):
+    """Calculates initial projection.
+
+    Calculates a projection at time step 0 and adds ghost cells on both
+    sides of the array.
+
+    Parameters
+    ----------
+    init_cond : InitialCondition object
+        Initial condition used for calculation.
+    mesh : Mesh object
+        Mesh for calculation.
+    basis: Basis object
+        Basis used for calculation.
+    quadrature: Quadrature object
+        Quadrature used for evaluation.
+    x_shift: float, optional
+        Shift in x-direction for each evaluation point. Default: 0.
+
+    Returns
+    -------
+    ndarray
+        Matrix containing projection of size (N+2, p+1) with N being the
+        number of mesh cells and p being the polynomial degree of the basis.
+
+    """
+    # Initialize output matrix
+    output_matrix = np.zeros((basis.polynomial_degree+1,
+                              mesh.num_cells+2*mesh.num_ghost_cells))
+
+    # 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[:, mesh.num_ghost_cells:
+                  output_matrix.shape[-1]-mesh.num_ghost_cells] = \
+        (basis_matrix * quadrature.weights) @ init_matrix
+
+    # Set ghost cells
+    output_matrix[:, :mesh.num_ghost_cells] = \
+        output_matrix[:, -2*mesh.num_ghost_cells:-mesh.num_ghost_cells]
+    output_matrix[:, output_matrix.shape[-1]-mesh.num_ghost_cells:] = \
+        output_matrix[:, mesh.num_ghost_cells:2*mesh.num_ghost_cells]
+
+    return output_matrix