Doesn't that our test functions should have enough differentiability. E.x., if the solution of PDE is C^2, then the test functions should be C^1? Thanks for help!

Hao, the discontinuous Galerkin method allows you to achieve approximation of $C^2$ functions by basis functions which are not globally continuous or differentiable. You write down the variational formulation of your PDE, taking care to account for jumps in the test functions and their derivatives over each of the geometric elements in your domain. Your discrete bilinear form thus includes volumetric integrals, as well as integrals along the interelement boundaries. You thus have to solve for not only the coefficients of the basis functions, but also their jumps. If you've set things up correctly, if your solution is globally $C^2$ then the jumps should all be zero, rendering your numerical approximation also $C^2$. There's considerable research into what formulations are stable, and in which contexts. You should look at papers by, for example, Doug Arnold, Bernardo Cockburn, and Dominik Schoetzau (just to mention a few). If you want a quick overview, http://www.cfm.brown.edu/people/jansh/resources/Publications/Lectures/RMMC08I.pdf The power of the DG method lies in its ability to achieve approximations for {\it nonsmooth} objects, but it's fine (if inefficient?) to use for elliptic PDE on smooth domains as well. 

