Why are identical test and basis functions preferred when using DG In discontinuous Galerkin FEM (as in continuous FEM) you begin by deriving the weak form of the PDE you are attempting to solve. The weak formulation guarantees exact conservation of the function you are trying to solve multiplied by a set of linearly independent functions (termed test functions). 
The next step is to expand the solution you are solving for as a linear combination of another set of linearly independent functions (termed basis functions). 
In typical DG method, the test and basis functions are the same i.e., they are derived from the same space. The only exception to this I have found are for PDEs with odd functions for which the same test and basis functions cannot be used.
Therefore, my question is why do people prefer using the same test and basis functions? Aside from it being the "traditional" method is there a mathematical reasoning? 
 A: Short answer: This is the way how the Galerkin procedure works: Test and basis functions are from the same function space (share the same basis).
Reason being: At some point in your weak formulation (assuming you apply DG to conservation laws) you have a term like
$$ \begin{align} \int_\Omega \boldsymbol{u}(x) \cdot \boldsymbol{v}(x) \mathrm d x =& \sum_i^N \int_{T_i} \sum_k^p \boldsymbol{u}_{i,k} \phi_{i,k}(x) \cdot  \sum_j^p \boldsymbol{e}_l \phi_{i,j}(x) \mathrm d x & \forall l \in 1, \dots, m\\
=& \sum_i^N  \int_{T_i} \sum_{k,j}^p \boldsymbol{u}_{i,k} \cdot \boldsymbol{e}_l \phi_{i,k}(x) \phi_{i,j}(x) \mathrm d x & \forall l \in 1, \dots, m \\ 
=& \sum_i^N \sum_{k,j}^p \boldsymbol{u}_{i,k} \cdot \boldsymbol{e}_l \int_{T_i} \phi_{i,k}(x) \phi_{i,j}(x) \mathrm d x & \forall l \in 1, \dots, m \\
=& \sum_i^N \sum_{k,j}^p \boldsymbol{u}_{i,k} \cdot \boldsymbol{e}_l \delta_{jk} = \sum_i^N \sum_{k}^p \boldsymbol{u}_{i,k} \cdot  \boldsymbol{e}_l & \forall l \in 1, \dots, m \\
=& \sum_i^N \sum_{k}^p \sum_j^m u_{i,k,j} \delta_{lj} = \sum_i^N \sum_{k}^p u_{i,k,l} & \forall l \in 1, \dots, m \end{align}$$
Above, $N$ is the number of cells, $p$ the number of basis functions adn $m$ the dimensionality of $\boldsymbol{u}$.
i.e., the integrals of the weak formulations boil down to Kronecker Deltas. Same holds if you have a bilinearform like $ \int_\Omega \nabla u(x) \cdot \nabla v(x) \mathrm d x$ as occuring for elliptic/parabolic PDEs.
You can read more about this for instance here, or here or here ... the list goes on, there is plenty literature out there. (If you drop for the moment the "discontinuous" buzzword).
