Sum Factorization in FEM

Shape functions, $\varphi ,$$\phi,$ are the building blocks of the finite element method. They describe how to interpolate a nodal solution field in terms of the nodal values on an element, ${\mathbf{u}}^e$$\bold{u}^e$

and evaluating these shape functions is part of every finite element calculation. For quadrilaterals, the shape functions for two commonly-used elements are:

$\begin{array}{c}\varphi (\mathit{\xi }):=\left[\begin{array}{c}(1-\xi )(1-\eta )\\ (1+\xi )(1-\eta )\\ (1-\xi )(1+\eta )\\ (1+\xi )(1+\eta )\end{array}\right]\end{array}$$\qquad \qquad \phi(\boldsymbol{\xi}) \coloneqq \begin{bmatrix} (1 - \xi)(1 - \eta) \\ (1 + \xi)(1 - \eta) \\ (1 - \xi)(1 + \eta) \\ (1 + \xi)(1 + \eta) \end{bmatrix}$

$\begin{array}{c}\varphi (\mathit{\xi }):=\left[\begin{array}{c}\frac{1}{4}(1-\xi )\xi (1-\eta )\eta \\ \frac{1}{2}(1-{\xi }^2)(1-\eta )\eta \\ \frac{1}{4}(1+\xi )\xi (1-\eta )\eta \\ \frac{1}{2}(1-\xi )\xi (1-{\eta }^2)\\ (1-{\xi }^2)(1-{\eta }^2)\\ \frac{1}{2}(1+\xi )\xi (1-{\eta }^2)\\ \frac{1}{4}(1-\xi )\xi (1+\eta )\eta \\ \frac{1}{2}(1-{\xi }^2)(1+\eta )\eta \\ \frac{1}{4}(1+\xi )\xi (1+\eta )\eta \end{array}\right]\end{array}$$\qquad \phi(\boldsymbol{\xi}) \coloneqq \begin{bmatrix} \frac{1}{4} (1 - \xi) \, \xi \, (1 - \eta) \, \eta \\\frac{1}{2} (1 - \xi^2) (1 - \eta) \, \eta \\ \frac{1}{4} (1 + \xi) \, \xi \, (1 - \eta) \, \eta \\ \frac{1}{2} (1 - \xi) \, \xi \, (1 - \eta^2) \\ (1 - \xi^2) \, (1 - \eta^2) \\ \frac{1}{2} (1 + \xi) \, \xi \, (1 - \eta^2) \\ \frac{1}{4} (1 - \xi) \, \xi \, (1 + \eta) \, \eta \\ \frac{1}{2} (1 - \xi^2) (1 + \eta) \, \eta \\ \frac{1}{4} (1 + \xi) \, \xi \, (1 + \eta) \, \eta \\ \end{bmatrix}$

If we look closely, we can see that these shape functions can be factored into products of simpler 1D interpolating polynomials in each direction:

Or, equivalently, the shape function associated with node $(i,j)$$(i,j)$ in the element is given by:

So, if we want to do a 2-dimensional interpolation, we can express it in terms of the composition of two 1-dimensional interpolations:

or in 3D,

If it is also the case that our quadrature points can also be written as cartesian products of some 1D quadrature rule (e.g. ${\mathit{\xi }}^q=({\xi }_s,{\xi }_t)$$\boldsymbol{\xi}^q = (\xi_s, \xi_t)$ or ${\mathit{\xi }}^q=({\xi }_r,{\xi }_s,{\xi }_t)$$\boldsymbol{\xi}^q = (\xi_r, \xi_s, \xi_t)$), then applying this sum factorization idea is beneficial for two reasons:

1. it lets us calculate all the quadrature point values at the same time in an efficient way

   • in 2D, the naive algorithm of individually interpolating each quadrature point of a $q\times q$$q \times q$ rule from the $n\times n$$n \times n$ nodal values involves $O(q^2{n}^2)$$O(q^2 n^2)$ operations, while the sum-factorization approach requires only $O(nq(n+q))$$O(n q (n + q))$.

   • in 3D, the naive algorithm of individually interpolating each quadrature point of a $q\times q\times q$$q \times q \times q$ rule from the $n\times n\times n$$n \times n \times n$ nodal values takes $O(q^3{n}^3)$$O(q^3 n^3)$ operations, while the sum-factorization approach requires only $O(q{n}^3+q^2{n}^2+q^{3}n)$$O(q n^3 + q^2 n^2 + q^3 n)$.

2. we only need to store evaluations of $\psi$$\psi$ over the 1D quadrature rule

   • in 2D, there are $n^2$$n^2$ shape functions and $q^2$$q^2$ quadrature points, so the naive storage approach requires $n^2{q}^2$$n^2 q^2$ values.

   • in 3D, there are $n^3$$n^3$ shape functions and $q^3$$q^3$ quadrature points, so the naive storage approach requires $n^3{q}^3$$n^3 q^3$ values.

   • If we only need the 1D shape functions evaluated over the 1D quadrature points, we only have to store $nq$$n q$ values, regardless of dimension.

Summary

Together, these two benefits make sum factorization a worthwhile optimization. However it is important to note that it is "fragile" and can only be used when all of the following requirements are met:

• quadrilateral or hexahedral elements only

  • the basis functions for simplex elements are not compatible

• the nodes in those elements must be arranged as cartesian products of some 1D nodes

  • not all formulations on quads/hexes are compatible (e.g. serendipity)

• the quadrature points must be expressible as the cartesian product of a 1D quadrature points

  • many applications that use irregular quadrature rules (e.g. mortar methods, contact, cut cells) do not benefit

• the polynomial order of the elements needs to be "big" to be worth it

  • sum factorization is technically still possible for low order elements, but it isn't faster in practice

• the number of quadrature points needs to be "big" to be worth it

  • many applications use single-point quadrature or "under-integrate" certain terms, so they see no benefit