Efficient assembly of high order continuous and discontinuous finite element operators

There is a growing consensus that state of the art Finite Element/Finite Volume technology is and will remain too computationally expensive to achieve the necessary resolution, even at the rate that computational power increases. The requirement for high fidelity computations naturally leads us to consider methods which have a higher order of grid convergence than the classical (formal) second order provided by most industrial grade codes. This indicates that higher-order discretization methods will at some point replace the solvers of today, at least for part of their applications. Although the potential of high order methods has been clearly demonstrated in literature, their inefficiency has been stressed by many authors. Only few publications are really dedicated to the matter; most of these focus on efficient, if not simplified, quadrature (Atkins and Shu, 1998). It is our opinion that this is not a desirable solution for non-linear problems as the potential savings are not tremendous. And these methods usually degrade or limit accuracy which may lead to decoupling for non-linear weakly hyperbolic equations. Moreover, it is shown in many contributions that the accuracy of a high order method strongly depends on the accuracy of the geometrical discretization (Bassi and Rebay, 1997 ; Bernard et al., 2008). In other words, high order methods require high order meshes and full precision quadrature. The finite element analysis process is classically decomposed into two parts: (i) the assembly process and (ii) the resolution process. At high orders, it is easy IMUM-2010, MIT August 17-20, 2010to show that the number of operations for the evaluation of the element matrices grows quickly with p, the polynomial order used for the approximation (O(p 9 ) in 3D). There is indeed a point, i.e. at some high order of approximation, where the assembly process inevitably becomes more expensive than the resolution process. In this work, we will show that it is possible to dramatically enhance the computational efficiency of the evaluation of finite element matrices by re-casting most of the floating point operations as large matrix-matrix multiplications. In this assembly process, no approximations are made on the quadrature or on the shape of the elements. The method can be applied both to continuous and discontinuous Galerkin formulations of systems of nonlinear PDEs. References H. L. Atkins, and C.-W. Shu, “Quadrature-free implementation of discontinuous Galerkin methods for hypebolic equations”, AIAA Journal , v. 35, p. 775-782, 1998. F. Bassi, and S. Rebay, “A high-order accurate discontinuous finite element solution of the 2D Euler equations”, Journal of Computational Physics, v. 138, p. 251-285, 1997. P.-E. Bernard, J.-F. Remacle, and V. Legat, “Boundary discretization for high order discontinuous Galerkin computations of tidal flows around shallow water islands”, International Journal for Numerical Methods in Fluids, v. 59, p. 535- 557, 2008