Matrices in two dimensions

[mrhardman]

This PR contributes two related features.

• The ability to calculate matrices for weak forms where two-dimensional integrals do not separate into a product of one-dimensional integrals. This occurs where the kernel $\rho(x_1,x_2)$ is not a separable function. Noting that $P_i(x_1)$, $R_i(x_2)$, $Q_i(x_1)$, $S_i(x_2)$ are from the set of Lagrange polynomials basis
  functions (and their first derivatives) in $x_1$ and $x_2$, an example of a non-separable matrix is

$$K_{ikjn} = \int^{x_{2u}}_{x_{2l}}\int^{x_{1u}}_{x_{1l}} P_i(x_1) R_k(x_2) Q_j(x_1) S_n(x_2)\rho(x_1,x_2) d x_1 d x_2.$$

• The syntax to create is a particular choice of $K_{ikjn}$=K[i,k,j,n] (noting the order of the indices) is

l_i = lagrange_x
l_j = lagrange_x
l_k = d_lagrange_dx
l_n = d_lagrange_dx
K_2D = finite_element_matrix(l_i,l_j,coordinate_x1,
                            l_k,l_n,coordinate_x2;
                            kernel_function=((x1,x2)->rho(x1,x2)),
                            max_iterations=10, quadrature_increment=10,
                            atol=1.0e-13, rtol=1.0e-13)

• A simple implementation of p-adaptive Gaussian quadrature using repeated calls to an 'internal' function _finite_element_matrix(), which calls FastGaussQuadrature to create the quadrature rule for a given number of quadrature points. This feature is necessary to make the caculation of $K_{ikjn}$ reliable for arbitrary $\rho(x_1,x_2)$.
• Additional tests.
• Updates to the README.md.