Consider the finite element method for the approximate solution of the equation $$ \nabla\cdot(a(x)\nabla u(x))=f(x) \text{ (in some domain)} $$ or, in weak form (assuming the right boundary conditions): Find $u\in H^1$ such that $$ \int a(x) \nabla u(x)\nabla v(x)=\int f(x) v(x)\quad\forall v\in H^1 $$
The finite element method replaces $u\in H^1$ and $v\in H^1$ by $u\in V$ and $v\in V$ for some finite dimensional vector space $V$ with a basis $(\psi_i)_{i=1}^n$ of compactly supported piecewise polynomials. This reduces the problem to a finite-dimensional linear system $$ Ax=b $$ where the entries of $A$ are $a_{ij}=\int a(x)\nabla \psi_i(x)\nabla\psi_j(x)$. When $a(x)$ is constant, these integrals can usually be computed exactly.
My question: When $a(x)$ is nonconstant, but $a\in C^m$, how exact do I need to approximate the entries $a_{ij}$? For example, if I use piecewise linear finite elements, and I want to maintain the convergence rate $\|u_h-u\|_{L^2}\leq Ch^2$ that exact computation of the $a_{ij}$ would provide, does it suffice to approximate the integrals by a midpoint rule?
There is a result in the Ciarlet book (Thm 4.1.6), which says that if the space contains polynomials of degree $k$, the numerical integration is exact for polynomials of degree $2k-2$, then we have $$ \|\nabla u-\nabla u_h\|_{L^2(\Omega)} \le c \ h^k \ ( \|a\|_{W^{k,\infty}}\|u\|_{H^{k+1}(\Omega)} + \|f\|_{W^{k,q}(\Omega)}) $$ where $q>2$, $k>n/q$.
Hence, for piecewise linear polynomials the midpoint rule is enough. I am not sure, whether this is sufficient to get the optimal rate for the $L^2$-error.