Implement interior penalty DG

Use this bilinear form:

Add a new file `grudge/models/elliptic.py` that realizes the above bilinear form. Allow for inhomogeneous Dirichlet and Neumann boundary conditions. Test for $h$-convergence to a known solution (quadratic polynomial on a square) for both boundary conditions.