Differentiation

A function represented at the quadrature points can be considered as a sum: if f_i,j is the (i,j)^th physical value of the function then the interpolating polynomial is:

$$\tilde{f}(a,b) = \sum_i^N \sum_j^N f_{ij} h^a_i(a) h^b_j(b)$$

Here h^a_i is the i^th interpolating polynomial for the quadrature chosen in the `a' direction. Thus, we can evaluate the derivative of this interpolating polynomial in each of the (a,b) directions.

$$\partial_a\tilde{f}(a,b) = \ \sum_n^N \sum_m^N f_{n,m} \partial_a(h^a_n)(a) h^b_m(b)$$

We can now evaluate the (i,j)^th quadrature value of the derivative as a matrix vector product:

$$\partial_a\tilde{f}(a_i,b_j) = \ \sum_n^N f_{nj} \partial_a(h^a_n)(a_i)$$

Repeating the above we can evaluate the `b' derivative (and the `c' derivative in three-dimensions). To calculate the x,y or z derivatives we simply apply the chain rule. In implementation we store the geometric factors for all the elements. This is more memory efficient for straight-sided triangles and tetrahedra since the geometric factors are a set of constants (except for (1-b) and (1-b)(1-c)² factors).