Which way is the right way to compute the integrals in finite element methods?

Question

Finite element methods involve integrals of functions that are not polynomials, and these integrals must be computed numerically.

For example, suppose that $f$ is the right-hand side of a Poisson problem and we use Lagrange elements of degree $k$. The right-hand side of the discrete linear system is computed from integrals

$\int_T f(x) p(x)$

where $p(x)$ is a polynomial of degree at most $k$ over the element $T$. What is the right way of computing that integral?

1. Use numerical quadrature for the integral. Considering the case $f = constant$ suggests that the quadrature formula should exactness at least $k$.

2. Replace the function $f(x)p(x)$ by a polynomial interpolant $\Pi[fp]$ and compute $\int_T \Pi[fp](x)$ exactly. Again, the interpolant should have at least degree $k$.

3. Replace the function $f(x)$ by a polynomial interpolant $\Pi[f](x)$ and compute $\int_T \Pi[f](x) p(x)$ exactly.

Both quadrature rules and interpolation require point evaluations of $f$. However, finite element textbooks always review quadrature rules. What are the considerations here?

Answer 1

Quadrature and replacing the integrand by a polynomial are identical. That is how quadrature rules are derived.

To see why this is true remember that quadrature evaluates the integrand at only a select number of points. When you do quadrature, you cannot distinguish between function $f(x)p(x)$ and $\Pi[fp](x)$ (where $\Pi[fp]$ is the function that interpolates $fp$ and $\Pi[fp](x)$ is that function evaluated at $x$) in the sense that $$ \int_a^b f(x) p(x) \, dx = \int_a^b \Pi[fp](x) \, dx $$ if you approximate both sides by quadrature. That's because during quadrature, you evaluate the integrands only at the interpolation points, and at these points the two integrands are the same by definition.

As a consequence, your approaches 1 and 2 are actually the same. Approach 3 is the same as the first two if you choose your quadrature points and weights appropriately.

Answer 2

Personally I typically use 1, and something similar to 2/3 though it'd be more correct to say that I decompose $f$ into a polynomial series of degree $k$ and then integrate that multiplied by $p$ exactly. I think there might be some way to prove some sort of equivalence with the interpolation methods you described.

I tend to find that 1 gives me the most flexibility and for a nodal FEM scheme where my nodes are located at the quadrature points your method will still converge at the expected rates. If I really need to I can interpolate to a higher order quadrature scheme's nodes and use the same fundamental implementation with higher order quadrature. This is how I tend to measure L-2 convergence since in that case you need a quadrature method of at least $O(2k+1)$, and I typically go a little bit higher just in case.

The other method I use of decomposing $f$ into a basis of degree $k$ tends to work well when doing something like an explicit time integration scheme because then you can pre-compute your inverse mass matrix and multiply it with the other matrix of coefficients to get a simple matrix-vector product to solve for your fields: $$ \partial_t q + \partial_x f = 0\\ \int_T p \partial_t q + \int_{\partial T} (f n)^* p - \int_T f \partial_x p = 0\\ \partial_t q_h = M^{-1} D f_h - M^{-1} L (f n)_h^*\\ M_{ij} = \int_T p_i(x) v_j(x)\\ D_{ij} = \int_T p_i(x) \partial_x v_j(x)\\ L_{ij} = \int_{\partial T} p_i(x) v_j(x) $$ You can of course do something similar with a second-order spatial derivative.

If I want to measure L-2 convergence I still need to perform higher-order quadrature, though.