A Large-Scale Comparison of Tetrahedral and Hexahedral Elements for Solving Elliptic PDEs with the Finite Element Method

TL;DR

This paper systematically compares tetrahedral and hexahedral finite elements across a large set of benchmark problems involving elliptic PDEs, highlighting their practical performance differences on automatically meshed 3D geometries.

Contribution

It introduces a comprehensive benchmark suite for evaluating FEM element types on real-world geometries, providing insights into their practical performance and system-level impacts.

Findings

Hexahedral elements often outperform tetrahedral elements in certain PDE solutions.

Mesh quality significantly influences FEM performance and accuracy.

The benchmark enables systematic comparison of FEM pipelines from meshing to solving.

Abstract

The Finite Element Method (FEM) is widely used to solve discrete Partial Differential Equations (PDEs) in engineering and graphics applications. The popularity of FEM led to the development of a large family of variants, most of which require a tetrahedral or hexahedral mesh to construct the basis. While the theoretical properties of FEM basis (such as convergence rate, stability, etc.) are well understood under specific assumptions on the mesh quality, their practical performance, influenced both by the choice of the basis construction and quality of mesh generation, have not been systematically documented for large collections of automatically meshed 3D geometries. We introduce a set of benchmark problems involving most commonly solved elliptic PDEs, starting from simple cases with an analytical solution, moving to commonly used test problem setups, and using manufactured solutions for thousands of real-world, automatically meshed geometries. For all these cases, we use state-of-the-art meshing tools to create both tetrahedral and hexahedral meshes, and compare the performance of different element types for common elliptic PDEs. The goal of his benchmark is to enable comparison of complete FEM pipelines, from mesh generation to algebraic solver, and exploration of relative impact of different factors on the overall system performance.

Tables

Table 1. Table 1. Comparison of performance of tetrahedral and hexahedral elements on several measures: time, memory, DOF and error, with one of the measures matched (marked in gray): for the first row of each comparison, we match time, second memory, etc. The best-performing (according to each measure) element is shown in green. For instance, by comparing P 2 subscript 𝑃 2 P_{2} with Q 1 subscript 𝑄 1 Q_{1} for the same error (last section of the table), P 2 subscript 𝑃 2 P_{2} is faster (fist column), it uses less memory (second column), and it has less DOFs (third column).

		time (s)	memory (MB)	DOF	error
$P_1$/$Q_1$	time	1.01 / 0.98	125 / 132	8,258 / 6,413	1.98e-03 / 6.93e-04
	memory	1.01 / 0.91	125 / 125	8,258 / 6,050	1.98e-03 / 7.49e-04
	DOF	1.01 / 1.07	125 / 149	8,258 / 7,139	1.98e-03 / 6.06e-04
	error	1.01 / 0.12	125 / 18	8,258 / 1,224	1.98e-03 / 1.87e-03
$P_2$/$Q_2$	time	16.86 / 17.56	2,236 / 2,033	59,885 / 44,541	1.85e-05 / 1.24e-05
	memory	16.86 / 18.77	2,236 / 2,241	59,885 / 48,951	1.85e-05 / 8.40e-06
	DOF	16.86 / 24.44	2,236 / 2,988	59,885 / 59,777	1.85e-05 / 5.43e-06
	error	16.86 / 11.22	2,236 / 1,451	59,885 / 35,017	1.85e-05 / 1.70e-05
$P_2$/$Q_1$	time	16.86 / 16.74	2,236 / 2,630	59,885 / 58,719	1.85e-05 / 1.58e-04
	memory	16.86 / 14.26	2,236 / 2,226	59,885 / 52,272	1.85e-05 / 1.70e-04
	DOF	16.86 / 17.52	2,236 / 2,669	59,885 / 59,777	1.85e-05 / 1.54e-04
	error	16.86 / 170.29	2,236 / 11,805	59,885 / 180,774	1.85e-05 / 6.36e-05

Equations

$$\mathcal{F}(x,u,\nabla u,D^{2}u)=b,\text{ subject to}$$

$$u=d\text{ on }\partial {\Omega }_D\text{and}\nabla u\cdot n=f\text{ on }\partial {\Omega }_N$$

$$u:\Omega \to {\mathbb{R}}^n,$$

$$\{\begin{array}{ll}-\Delta u=b& \text{on }\Omega \\ u=d& \text{on }\partial {\Omega }_D\\ \nabla u\cdot n=f& \text{on }\partial {\Omega }_N.\end{array}$$

$$\{\begin{array}{ll}-\mu \Delta u+\nabla p=b& \text{on }\Omega \\ -\text{divergence}u=0& \text{on }\Omega \\ u=d& \text{on }\partial {\Omega }_D\\ \mu (\nabla u+{\nabla }^{T}u)\cdot n-pn=f& \text{on }\partial {\Omega }_N\end{array}$$

$$\{\begin{array}{ll}-\text{divergence}\sigma [u]=b& \text{on }\Omega \\ u=d& \text{on }\partial {\Omega }_D\\ \sigma [u]n=f& \text{on }\partial {\Omega }_N.\end{array}$$

$${\sigma }^L[u]=2\mu ϵ[u]+\lambda \text{tr}ϵ[u]Iϵ[u]=\frac{1}{2}\left(\nabla u^T+\nabla u\right),$$

$$\mu =\frac{E}{2(1+\nu )},{\lambda }_{3D}=\frac{E\nu }{(1+\nu )(1-2\nu )},\mathrm{and}{\lambda }_{2D}=\frac{\nu E}{1-{\nu }^2}.$$

$$\{\begin{array}{ll}-\text{divergence}(2\mu ϵ[u]+pI)=b& \text{on }\Omega \\ \text{divergence}u-{\lambda }^{-1}p=0& \text{on }\Omega \\ u=d& \text{on }\partial {\Omega }_D\\ {\sigma }^N[u]\cdot n=f& \text{on }\partial {\Omega }_N\end{array}$$

$$\sigma [u]=\mu (F[u]-F[u{]}^{-T})+\lambda \ln (\det F[u])F[u{]}^{-T}F[u]=\nabla u+I,$$

$$\begin{array}{rl}{S}_{\text{2D}}^2=& {\sigma }_{0,0}^2-{\sigma }_{0,0}{\sigma }_{1,1}+{\sigma }_{1,1}^2+3{\sigma }_{0,1}{\sigma }_{1,0}\\ S_{\text{3D}}^2=& \frac{({\sigma }_{0,0}-{\sigma }_{1,1}{)}^2+({\sigma }_{2,2}-{\sigma }_{1,1}{)}^2+({\sigma }_{2,2}-{\sigma }_{0,0}{)}^2}{2}+\\ & 3({\sigma }_{0,1}{\sigma }_{1,0}+{\sigma }_{2,1}{\sigma }_{1,2}+{\sigma }_{2,0}{\sigma }_{0,2}).\end{array}$$

$$u$$

$$3/4e^{-((9x_1-2{)}^2+(9x_2-2{)}^2+(9x_3-2{)}^2)/4}$$

$$+3/4e^{-(9x_1+1{)}^2/49-(9x_2+1)/10-(9x_3+1)/10}$$

$$+1/2e^{-((9x_1-7{)}^2+(9x_2-3{)}^2+(9x_3-5{)}^2)/4}$$

$$-1/5e^{-(9x_1-4{)}^2-(9x_2-7{)}^2-(9x_3-5{)}^2},$$

$$u(x_1,x_2,x_3)=\frac{1}{80}\left(\begin{array}{c}{x}_1{x}_2+x_1^2+x_2^3+6x_3\\ x_1{x}_3-x_3^3+x_1{x}_2^2+3x_1^4\\ x_1{x}_2{x}_3+x_2^2{x}_3^2-2x_1,\end{array}\right)$$

Full text

We introduce a set of benchmark problems involving most commonly solved elliptic PDEs, starting from simple cases with an analytical solution, moving to commonly used test problem setups, and using manufactured solutions for thousands of real-world, automatically meshed geometries. For all these cases, we use state-of-the-art meshing tools to create both tetrahedral and hexahedral meshes, and compare the performance of different element types for common elliptic PDEs.

The goal of his benchmark is to enable comparison of complete FEM pipelines, from mesh generation to algebraic solver, and exploration of relative impact of different factors on the overall system performance.

As a specific application of our geometry and benchmark dataset, we explore the question of relative advantages of unstructured (triangular/tetrahedral) and structured (quadrilateral/hexahedral) discretizations. We observe that for Lagrange-type elements, while linear tetrahedral elements perform poorly, quadratic tetrahedral elements perform equally well or outperform hexahedral elements for our set of problems and currently available mesh generation algorithms. This observation suggests that for common problems in structural analysis, thermal analysis, and low Reynolds number flows, high-quality results can be obtained with unstructured tetrahedral meshes, which can be created robustly and automatically.

We release the description of the benchmark problems, meshes, and reference implementation of our testing infrastructure to enable statistically significant comparisons between different FE methods, which we hope will be helpful in the development of new meshing and FEA techniques.

Tetrahedral Mesh, Hexahedral Mesh, Finite Element, Comparison

††copyright: acmcopyright††journal: TOG††journalyear: 2022††journalvolume: 1††journalnumber: 1††article: 1††publicationmonth: 1††price: 15.00††doi: 10.1145/3508372††ccs: Computing methodologies Mesh models††ccs: Computing methodologies Mesh geometry models††ccs: Computing methodologies Volumetric models††ccs: Computing methodologies Simulation evaluation

1. Introduction

The finite element method (FEM) is commonly used to discretize partial differential equations (PDEs), due to its generality, rich selection of elements adapted to specific problem types, and wide availability of commercial implementations. At a high level, a FE analysis code takes as input the domain boundary, the boundary conditions, and the governing equations of the phenomena of interest, and computes the solution everywhere in the domain.

As an initial step in this procedure, the domain typically has to be discretized in a finite collection of elements. Many choices are possible, ranging from unstructured grids of tetrahedra to perfectly regular grids of cubes. Despite the large amount of research on mesh generation, we were unable to find a systematic study answering a basic question: “What are the practical pros and cons of using unstructured (triangular/tetrahedral) or structured (quadrilateral/hexahedral/grids) discretizations for commonly used elliptic PDEs?”.

This question is critical to inform the development of meshing algorithms: while tetrahedral meshes are easier to generate automatically, hexahedral meshes (i.e., meshes that are composed of only deformed cubes) are much more difficult to adapt to objects with complex geometries, while maintaining high mesh quality. One of the arguments motivating development of these more complex algorithms is a common belief that hexahedral elements yield better accuracy for a given computational cost (see the introduction of, e.g., (Lyon et al., 2016; Guo et al., 2020; Bernard et al., 2016)).

The overall aim of our work is to provide an extensive benchmark for comparing the performance of FE pipelines, including automatic meshing, FE basis construction, and algebraic system solvers, on a set of most common elliptic PDEs and a set of realistic geometries. As an immediate application, we explore the performance of widely used families of elements, coupled with standard solvers, on a large set of meshes generated using currently available meshing algorithms.

More specifically, we compare the efficiency of different elements, that is, how much time is typically required to obtain a solution with a given accuracy for different element types on automatically generated unstructured meshes, on manually and automatically generated semi-structured meshes, and on regular lattices.

We consider standard Lagrangian bases (Szabó and Babuška, 1991; Ciarlet, 2002a) of varying degrees, as well as serendipity (Zienkiewicz et al., 2005) elements (for hexahedra only), which are by far the most popular brick element. Finally, we perform several comparisons using spline-based elements (Hughes et al., 2005), which have recently gained popularity in the IsoGeometric Analysis (IGA) community. While this clearly does not reflect the broad range of existing element types and PDEs in the literature, it includes the most popular general-purpose elements currently used in commercial and open-source FE systems. For solving the resulting linear systems, we consider both state of the art direct (De Coninck et al., 2016) and iterative solvers (Falgout and Yang, 2002).

We collected a set of test problems of varying complexity for elliptic PDEs (including Poisson, linear elasticity, Neo-Hookean elasticity, incompressible elasticity, and incompressible Stokes equations). Our set includes common simple test problems (where most of the hex-meshes are grids): beam bending, beam twisting, driven cavity flow, planar domain with a hole, elasticity problems with singular solutions, as well as a large-scale benchmark of manufactured solutions (Salari and Knupp, 2000) on $3\,200$ automatically meshed, real-world, complex 3D models. Our model collection includes both CAD models and scanned geometries, providing a realistic sampling of analysis scenarios. We use TetWild (Hu et al., 2018; Hu et al., 2020) and MeshGems (Spatial, 2018) to generate the tetrahedral and hexahedral meshes respectively (we also included the state-of-the-art meshes from Hexalab (Bracci et al., 2019)).

This combination of test models, 3D meshes for these models, elements and solvers is representative of many common FE application scenarios.

We quantify (to the best of our knowledge, for the first time) the overall performance differences between these two families of elements. Our main conclusion is that, while linear elements on triangular/tetrahedral meshes exhibit well-known problems, quadratic tetrahedral elements perform similarly or better (i.e., require similar or less time to compute a solution with a given accuracy) than Lagrangian elements on semi-structured hexahedral meshes, and are somewhat inferior (but still competitive, especially considering tetrahedral meshing is much faster and more robust) to the performance of spline elements on regular lattices when a direct solver is used. Combined with available state-of-the-art robust meshing techniques, quadratic tetrahedral elements are a good choice to realize a fully automatic pipeline,.e.g., for SciML applications, or shape optimization, without sacrificing performance compared to hexahedral elements, which require far more complex and less robust mesh generation. More detailed conclusions are presented in Section 6.

We emphasize that our study is limited to a specific set of PDEs, commonly used geometry-agnostic linear solvers, and state-of-the-art meshing algorithms; we leave adding dynamic scenarios, multi-physics, different linear solvers, and other extensions as future work – the provided framework can be readily extended to these cases. We also note that adaptive refinement is simpler for hexahedral meshes and, as a consequence, adaptive geometric multigrid solvers are more readily available (Alzetta et al., 2018), although it is possible to develop similar solvers for tetrahedral meshes (Kohl et al., 2019). While the outcome of our study should not be interpreted as a reason to favor tetrahedral discretizations in all situations (and there are applications of hexahedral meshes outside the scope of FEM discretizations, such as lattice structure design), it does point to the need for direct experimental evaluation of meshing strategies, in the context of specific target applications.

We provide the complete source code111https://github.com/polyfem/polyfem/ for the integrated analysis pipelines we tested, the dataset we used222https://archive.nyu.edu/handle/2451/44221, the benchmark solutions, and the scripts to reproduce all results333https://github.com/polyfem/tet-vs-hex, to enable researchers and practitioners to easily expand this study to additional mesh types (such as polyhedral meshes) and bases.

This study is divided into five sections: we first introduce the closest related works on meshing and analysis (Section 2). We then overview the background on mesh types, basis, and the model PDE that we consider in this study (Section 3). We divide the experimental evaluation into a set of individual experiments, targeting a set of common test problems (including problems with singularities) in Section 4, and then perform a large scale analysis on thousands of automatically generated meshes in Section 5. We finally draw conclusions and identify open challenges in Section 6.

2. Related Work

We first review existing comparisons of different types of finite elements (Section 2.1), then briefly discuss commonly used finite element software (Section 2.2) and the state-of-the-art meshing algorithms (Section 2.3).

2.1. FEA on Unstructured and Structured Meshes

To the best of our knowledge, our study is the first large-scale comparison between different commonly used types of elements in FEM. However, there are multiple existing comparisons focused on specific models and physics.

In (Cifuentes and Kalbag, 1992), the authors conclude that quadratic tetrahedral meshes lead to roughly the same accuracy and time as linear hexahedral meshes, by comparing solutions for several simple structural problems. By evaluating the eigenvalues of the stiffness matrices of various nonlinear and elastoplastic problems, (Benzley et al., 1995) reports that, in their study, linear hexahedral meshes are superior to linear tetrahedral meshes. The authors also show that linear hexahedral meshes are slightly superior to quadratic tetrahedral meshes in the nonlinear elastoplastic analysis experiment.

A more recent work, (Tadepalli et al., 2010, 2011), focuses on modeling footwear with a nonlinear incompressible material model under shear force loading conditions. The conclusion of these works is that trilinear hexahedral meshes are superior to linear tetrahedral meshes, and that quadratic tetrahedral elements are computationally more expensive compared to trilinear hexahedral elements, but have higher accuracy. (Wang et al., 2004) compares tetrahedral and hexahedral meshes on linear static problems, modal and nonlinear analysis. The study concludes that quadratic tetrahedral and hexahedral elements have similar performance, but quadratic hexahedra are computationally more expensive. The same study also confirms that linear tetrahedra are too stiff for large deformations, and linear hexahedra with large corner angles should be avoided in regions of stress concentration. The study is restricted to a small set of geometries and focuses on manual hexahedral mesh generation. Our study instead focuses on automatic meshing algorithms for both tetrahedral and hexahedral meshes, and we provide experimental results on thousand of complex geometric models and a wide array of elliptic PDEs.

In medical applications, results for femur models (Ramos and Simões, 2006) show that linear tetrahedral meshes of the simplified femur model lead to a closer agreement with the theoretical ones, while quadratic hexahedral meshes are more stable and the result is less affected by mesh refinement. On a kidney model, (Bourdin et al., 2007) observes that both linear and quadratic tetrahedral meshes are slightly stiffer than hexahedral meshes, but are more stable when high impact energies are present in the simulation. For heart mechanics and electrophysiology, (Oliveira and Sundnes, 2016) notes that quadratic hexahedra are slightly better than quadratic tetrahedra in the mechanics regime, while linear tetrahedral meshes are the best choice for the electrophysiology problem.

2.2. Finite Element Analysis Software

There exists a large number of libraries and software for finite-element analysis, both open-source and commercial. An exhaustive comparison of all existing packages444A non-exhaustive list of open-source FEA packages known to the authors include, in alphabetical order, code_aster (EDF, 2018), Deal.II (Alzetta et al., 2018), DOLFIN (FEniCS) (Alnæs et al., 2015), ElmerFEM (Elmer, 2018), FEATool Multiphysics (MATLAB) (Ltd., 2019), Feel++ (Prud’homme et al., 2012), FEI (Trilinos) (Heroux et al., 2005), Firedrake (McRae et al., 2016), FreeFEM++ (Hecht, 2012), GetDP (Geuzaine, 2008), GetFEM++ (Renard and Pommier, 2018), libMESH (Kirk et al., 2006), MFEM (MFEM, 2020), Nektar++ (Cantwell et al., 2015), NGSolve (Schöberl, 2014), OOFEM (Patzák, 2012), PolyFEM (Schneider et al., 2019b), Range (Šoltys, 2019), SOFA (Faure et al., 2012), and VegaFEM (Barbič et al., 2012). is beyond the scope of this paper, therefore we discuss only several representative packages. We point out an interesting project (Ladutenko, 2018) attempting to maintain a complete list of FEA packages with a list of characteristics.

Our goal is to investigate and compare the performance of FEM on meshes with tetrahedral and hexahedral elements, using the standard Lagrangian basis functions and serendipity elements commonly used in engineering applications, as well as spline elements used in IGA.

Open-source packages such as FEniCS (Alnæs et al., 2015), GetFEM++ (Renard and Pommier, 2018), libMesh (Kirk et al., 2006), and MFEM (MFEM, 2020) support both tetrahedral and hexahedral meshes, although very few (e.g., libMesh) implement both the 20-(serendipity) and 27-nodes variant for quadratic hexahedral elements. Deal.II (Alzetta et al., 2018) is another popular open-source FEA library, however it only supports quadrilateral and hexahedral elements. Commercials packages such as ANSYS (ANSYS Inc., 2019), Abaqus (ABAQUS Inc., 2019), COMSOL Multiphysics (COMSOL Inc., 2018) support Lagrangian tetrahedral elements, but surprisingly often implement only serendipity elements for hexahedra (Zienkiewicz et al., 2005, Chapter 6). Given their popularity, we included serendipity elements in our study in addition to traditional Lagrangian elements.

Another increasingly popular choice of bases for hexahedral meshes are B-splines and NURBS, most commonly used in the context of isogeometric analysis (IGA) (Hughes et al., 2005). The popularity of spline bases stems from the fact that they have only one dof per element independently of the degree (however, the support of each basis function grows accordingly, and, as a consequence the stiffness matrices become less sparse). Defining this type of element on fully general hexahedral domains is an open problem (Aigner et al., 2009; Martin and Cohen, 2010; Li et al., 2013). Due to their rising popularity, we deem important to include experiments with these elements in our study, but restrict them to cases where a regular lattice mesh is used.

Since none of these libraries implements both Lagrangian (tetrahedral and hexahedral), serendipity, and spline basis functions (hexahedral only) in the same framework, we added all the elements and basis used in this study to our own open-source FEA library (Schneider et al., 2019b) to ensure a fair comparison. PolyFEM (Schneider et al., 2019b) supports all these element types and interfaces with Hypre (Falgout and Yang, 2002) and PARDISO (De Coninck et al., 2016; Verbosio et al., 2017; Kourounis et al., 2018) for the solver and Eigen (Guennebaud et al., 2010) for linear algebra.

2.3. Meshing

Three-dimensional mesh generation has been thoroughly studied in multiple communities (Shewchuk, 2012; Carey, 1997; Owen, 1998; Tautges, 2001). For the sake of brevity, we restrict our review to the techniques generating pure tetrahedral or pure hexahedral meshes, which are the focus of our study, with an emphasis on methods implemented in readily available open-source or commercial libraries.

Tetrahedral Meshing

The most efficient, popular, and well-studied family of algorithms tackles the generation of meshes satisfying the Delaunay condition (Chew, 1993; Shewchuk, 1996, 1998; Ruppert, 1995; Shewchuk, 2012; Sheehy, 2012; Remacle, 2017; Du and Wang, 2003; Alliez et al., 2005; Tournois et al., 2009; Murphy et al., 2001; Cohen-Steiner et al., 2002; Chew, 1987; Si and Gärtner, 2005; Shewchuk, 2002; Si and Shewchuk, 2014; Si, 2015; Cheng et al., 2008; Boissonnat and Oudot, 2005; Jamin et al., 2015; Dey and Levine, 2008; Chen and Xu, 2004). These methods are robust if the input is a point cloud, but might fail if the boundary of a shape has to be preserved exactly (Hu et al., 2018; Hu et al., 2020).

To overcome these robustness limitations, alternative approaches are based on a background grid (Molino et al., 2003; Bronson et al., 2013; Labelle and Shewchuk, 2007; Doran et al., 2013; Bridson and Doran, 2014). The idea is to fill the bounding box of the 3D input surface with either a uniform grid or an adaptive octree, whose convex cells are trivial to tetrahedralize. These methods achieve high quality in the interior of the mesh (where the grid is regular), but introduce badly shaped elements near the boundary, which is often the region of interest in many practical simulations. On the other hand, front-advancing methods (Cuillière et al., 2013; Alauzet and Marcum, 2014; Haimes, 2014) start by marching from the boundary to the interior, adding one element at a time, pushing the problematic elements into the interior where the advancing fronts meet.

All these methods are unable to handle commonly occurring input surfaces which contain degenerated faces, gaps, and self-intersections. These types of defects are, unfortunately, common in CAD models, due to the NURBS representation (with a fixed degree) not being closed under boolean operations. To the best of our knowledge, the only method that was demonstrated to be capable of handling these cases robustly is TetWild (Hu et al., 2018). It is based on a hybrid numerical representation to ensure correctness, and it allows a small, controlled deviation from the input surface to achieve a good element quality. We used this technique to generate all unstructured tetrahedral meshes in this study.

Hexahedral Meshing

aims at filling the volume enclosed by an input surface with hexahedra. Hexahedra also need to have a good shape to ensure good solution approximation. The natural tensor-product structure of a hexahedron enables to define tensor-product bases, and, e.g., use spline-based elements, but dramatically increases the complexity of meshing algorithms. Semi-manual or interactive approaches are usually employed, such as sweeping and advancing front methods (Shepherd and Johnson, 2008; Gao et al., 2016; Livesu et al., 2016), which are used in commercial software such as (ANSYS Inc., 2019; Coreform, 2020).

By allowing lower element quality, one can design automatic approaches based on regular lattices (Schneiders and Bünten, 1995; Schneiders, 1996; Su et al., 2004; Zhang and Bajaj, 2006; Zhang et al., 2007) or on octrees (Schneiders et al., 1996; Zhang and Bajaj, 2006; Ito et al., 2009; Maréchal, 2009; Qian and Zhang, 2010; Ebeida et al., 2011; Zhang et al., 2013; Elsheikh and Elsheikh, 2014; Owen et al., 2017; Spatial, 2018).

Polycube methods (Gregson et al., 2011; Livesu et al., 2013; Huang et al., 2014; Fang et al., 2016; Fu et al., 2016; Li et al., 2013; Zhao et al., 2018) and field-aligned parameterization-based methods (Nieser et al., 2011; Huang et al., 2011; Li et al., 2012; Jiang et al., 2014; Solomon et al., 2017; Liu et al., 2018) aim at producing hexahedral meshes with as few irregular edges and vertices as possible, but designing robust algorithms of this type is still an open problem. Sample results from some of the previous methods have been recently collected into a single repository (Bracci et al., 2019), which we use in our study. We also generate a new dataset composed of 3200 hexahedral meshes using the commercial MeshGems-Hexa software (Spatial, 2018).

3. Background

3.1. FEM bases

There is a multitude of different definitions of bases for both tetrahedral (or triangular) and hexahedral (or quadrilateral) element shapes, with different elements tailored to specific types of problems (e.g., axisymmetric elements, shell elements, plasticity elements, etc.). In our comparison, we target the most common choices: we use the standard linear and quadratic Lagrange bases for tetrahedra, which we denote $P_{1}$ and $P_{2}$ respectively, and hexahedra, with $Q_{1}$ denoting linear tensor-product basis and $Q_{2}$ quadratic tensor-product basis (Szabó and Babuška, 1991; Ciarlet, 2002b). We also use the serendipity basis (Zienkiewicz et al., 2005), commonly used in commercial software, and spline basis (Hughes et al., 2005) for hexahedral elements. We use the standard Galerkin formulation (Szabó and Babuška, 1991; Ciarlet, 2002b) with Gaussian quadrature for all our experiments, avoiding non-standard quadrature.

3.2. Mesh and solution characterization

We use the number of vertices as a measure of the resolution of tetrahedral and hexahedral meshes, as the number of vertices is often used by the meshing algorithms as the “budget” that the meshing algorithms can use to create the best possible mesh, and the number of vertices is equal to the number of degrees of freedom in the case of linear (or tri-linear) elements.

In addition to this particular choice, we also investigate other metrics for a specific example (Table 1), and provide an interactive plot that allows one to compare our results using 24 different measures: solution error measured using $H^{1}$, $H^{1}$ semi-norm, $L_{2}$, $L^{\infty}$, $L^{\infty}$ of gradient, and $L^{8}$ norms; mesh average edge length, minimum edge length and number of vertices; the system matrix size and the number of non zero entries, the numbers of basis functions, dofs, elements, and pressure basis functions; timings for loading mesh data, building basis functions, computing the right-hand side, assembling the system matrix, solving the system, computing the errors, total time and time without right-hand side assembly.

3.3. Model PDEs

We selected the following set of representative elliptic problems: (1) Poisson; (2) incompressible stationary Stokes fluid flow equations; (3) elasticity with linear Hooke’s law as the constitutive equation; (4) Neo-Hookean elasticity (5) incompressible linear elasticity. We list the corresponding PDEs for completeness.

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$ be the domain with boundary $\partial\Omega$. We aim to solve

[TABLE]

for the function

[TABLE]

where $D^{2}$ is the matrix of second derivatives, $b$ is the right-hand side, $\partial\Omega_{D}\subset\partial\Omega$ is the part of the boundary with Dirichlet boundary conditions, and $\partial\Omega_{N}\subset\partial\Omega$ is the part of the boundary with Neumann boundary conditions. Since we consider second-order PDEs only, $\partial\Omega_{D}\cap\partial\Omega_{N}=\emptyset$. The form of $\mathcal{F}$ and the role of $u$ depends on the specific PDE.

We consider polygonal and polyhedral domains $\partial\Omega$ (possibly non-convex). The right-hand side $b$ in our test examples is analytic, the boundary $d$ is continuous and piecewise-smooth, and $f$ is piecewise smooth (but possibly with finite-jump discontinuities); under these assumptions, the weak solutions of the equations we consider are (at least) continuous, but the solution derivatives may be singular. We primarily focus on the error in the solution itself, rather than the derivative error, although consider the stress for some elasticity examples. We state the model problems in the strong form, but only the weak solutions exist for many of the test cases.

Poisson Equation

This problem is given by

[TABLE]

Incompressible Steady Stokes Equations

The Stokes equations provide the relationship between the velocity $u$ and the pressure $p$ for an incompressible fluid with viscosity $\mu$.

[TABLE]

Elasticity

Elasticity PDEs are formulated in terms of the stress tensor $\sigma[u]$ (which depends on the displacement $u$) as

[TABLE]

In this case the right-hand side $b$ plays the role of a body force, the Dirichlet boundary conditions are fixed displacement, and the Neumann ones are surface tractions.

Material models define how the stress $\sigma$ is related to the displacement field $u$. For the linear Hookean model,

[TABLE]

where $\epsilon[u]$ is the strain tensor, $\lambda$ is the first Lamé parameter, and $\mu$ is the shear modulus. There are two common assumptions reducing the elasticity problem to a 2D problem, plane stress and plane strain; in our experiments we are using plane stress. In this case, the elasticity equation has the same form but with different constants (Hughes, 2012):

[TABLE]

Incompressible materials form a separate class: in 3D, an isotropic material has Poisson ratio equal to 0.5, and the previous equation is not well-defined, as $\lambda$ becomes infinite. While isotropic materials in plane stress state cannot have this problem, as the isotropic Poisson ratio cannot exceed 0.5, anisotropic materials can have Poisson ratio 1 for in-plane deformations, and thus can be 2D-incompressible, which geometrically corresponds to the area of the cross-section of a material element preserved under deformations  (Lee and Lakes, 1997)). As a consequence, equations for 2D-incompressible materials in plane stress state are also of interest. Additionally, when $\lambda$ grows, the linear system arising from the discretization of the PDE becomes unstable. A common way to avoid such problem is to introduce a Lagrange-multiplier-like function in the form of the pressure $p$. This leads to a mixed formulation of elasticity similar to Stokes equations which is stable for large $\lambda$s, and reduces to incompressible elasticity for $\lambda^{-1}\rightarrow 0$.

[TABLE]

Finally, in the Neo-Hookean material model the stress is a nonlinear function of strain.

[TABLE]

where $F[u]$ is the deformation gradient.

For elasticity problems, we often use the von Mises stresses

[TABLE]

Note that the stresses are discontinuous since they depend on the gradient of the displacement which is only $C^{0}$ for our discretizations. To mitigate visual artefacts we average the stresses around vertices in our plots.

3.4. Linear Solvers

All FEM problems we consider require to solve a linear system, which, as the mesh size grows, dominates the running time. A vast amount of research has been invested in developing efficient and robust linear solvers. In our study we use two state-of-the art solvers: Pardiso (De Coninck et al., 2016) a direct solver using the Cholesky factorization, which we use for smaller problems, and Hypre (Falgout and Yang, 2002) an algebraic multigrid solver, which we use for larger problems. Direct solvers work particularly well in 2D, but scale poorly for 3D problems. We leave as future work a more detailed study on the effect of the linear solver on the solution time. The conclusions of this study hold for both types of solvers for our experimental setup and test problems.

4. Common Test Problems

We collected a number of standard test cases to cover different physical phenomena and different scenarios: fluid simulation (Section 4.1), linear elastic time dependent (Section 4.2), linear elastic bars (Section 4.3), linear orthotropic material models (Section 4.4), meshes with high aspect-ratio for linear elastic bars (Section 4.5), classical plane with hole with symmetric boundary conditions for compressible and nearly incompressible material (Section 4.6), nearly incompressible linear material (Section 4.7), nonlinear Neo-Hookean material (Section 4.8), and nonlinear Neo-Hookean material with high stresses (Section 4.9).

Most of the solution domains are chosen to simplify manual creation of hexahedral meshes: the simulations will be performed on an unstructured tetrahedral mesh and a nearly regular lattice with the same number of vertices. Experiments in Sections 4.2 to 4.7 are run on a MacBook Pro 3.1GHz Intel Core i7, 16GB of RAM, and 8 threads. Experiments in Sections 4.8 and 4.9 are run on a cluster node with 2 Xeon E5-2690v4 2.6GHz CPUs and 250GB memory, each with max 128GB of reserved memory and 8 threads. For all experiments, we use the PolyFEM library (Schneider et al., 2019b), which uses the Pardiso (De Coninck et al., 2016; Verbosio et al., 2017; Kourounis et al., 2018) direct solver, and Newton iterations for the nonlinear problems.

Note that, for completeness, we also validated PolyFEM on the example in Figure 3 for linear and quadratic tetrahedra and serendipity hexahedra on Hooke material against Abaqus. The results are identical up to numerical precision.

4.1. Incompressible Stokes

We use a planar square domain mesh with 4 229 vertices for the triangle mesh and 4 225 vertices for the regular grid. We simulate the Stokesian fluid (2) with viscosity $\mu=1$ in the standard “driven cavity” example: the fluid has zero boundary conditions on 3 of the 4 sides and a tangential velocity of 0.25 on the left side. Figure 1 shows the results for mixed linear (for the pressure) and quadratic (for the velocity) elements: the results are indistinguishable between hexahedral and tetrahedral elements.

4.2. Time-Dependent Linear Elasticity

We consider the dynamics of a suspended object under gravity: we fix the top part of a unit square with material parameters $E=200$ and $\nu=0.35$ and apply a constant body force of $20$ in the $y$ direction. We integrate the dynamic simulation for $t$ from [math] to $0.5$ with 40 time steps integrated with Newmark (Newmark, 1959). We mesh the domain at a coarse and fine resolution, both for triangles and for quads. Figure 2 shows the displacement in the $x$ direction of the bottom left corner for the 4 discretizations, using linear and quadratic elements.

4.3. Transversally Loaded Beam

In this experiment, we consider beams with different cross-sections (square, circular, and I-like) in the $xy$-plane of length $L$. The beam is fixed (i.e., zero Dirichlet conditions are applied) at the end ($z=L$), and different tangential forces $f=[0,-f_{y},0]^{T}$, $f_{y}\in[-0.1,-2]$, are applied at $z=0$, opposite to the fixed side. The rest of the boundary is left free and we do not apply any body force. For these experiments we use linear isotropic material model (4) with Young’s modulus $E=210\,000$ and Poisson’s ratio $\nu=0.3$. We study the displacement at the bottom corner of the moving end ($z=0$) in the $y$ direction and compare it with a dense solution to compute the error $e$ (note that the solution is singular only at $z=L$, far from the evaluation points). We report as $e_{f}$ the slope of the linear fit of the error as a function of the force magnitude. We also report the basis construction time $t_{b}$, assembly time $t_{a}$, solve time $t_{s}$, and total time $t$. Note that all the timings reported are averaged over $10$ different runs per force sample.

Square Cross-section

For running the simulation, we use a square cross-section of side $s=20$, length $L=100$ and mesh it with a tetrahedral mesh with $739$ vertices and a hexahedral mesh (regular grid) with $750$ vertices. Figure 3 shows the errors compared with the dense solution, where trilinear hexahedral elements outperform linear tetrahedral elements but the quadratic counterparts are indistinguishable. Timing-wise, the quadratic tetrahedra are slightly better.

We created a sequence of hexahedral and tetrahedral meshes with similar errors for a force $f=[0,-2,0]^{T}$. Figure 4 shows that for a given error, $P_{2}$ discretization is around four times faster than $Q_{1}$, and $\mathrm{SPLINE}_{2}$ have a slight advantage over $P_{2}$. Note that both $Q_{1}$ and $\mathrm{SPLINE}_{2}$ are constructed over a perfectly regular grid, while the $P_{2}$ elements are defined over an unstructured tetrahedral mesh.

Finally, we created a sequence of hexahedral meshes that matches total time, total memory, total number of degrees of freedom, and error of the tetrahedral mesh in Figure 3 for both linear and quadratic elements. Table 1 summarizes our findings: $P_{1}$ is significantly worse than $Q_{1}$ but the two quadratic discretizations produce similar results. $P_{2}$ overall performs better than $Q_{1}$.

Circular Cross-section

We consider a beam of length $L=100$ with a circular cross-section of diameter $d=20$. We created a tetrahedral mesh with $2\,252$ vertices and a hexahedral mesh with $2\,288$ vertices (note that in this case the mesh is not a regular grid anymore), by extruding a quad mesh generated with (Jakob et al., 2015). Figure 5 shows similar $y$-displacement errors as for the square cross-section, $P_{1}$ produces low-quality results, while $P_{2}$ and $Q_{2}$ are similar.

I-beam Cross-section

We use an I-beam (the bounding box of the cross-section is $125\times 154$) of length $L=473.11$. The tetrahedral mesh has $6\,102$ while the hexahedral mesh, generated by extruding a quad mesh, has $6\,080$ vertices, results are shown in Figure 6.

4.4. Orthotropic Material

We repeated the previous experiment using linear orthotropic material parameters (carbon fiber). The material parameters are obtained from (Pardini and Gregori, 2010). Three Young moduli are 167, 33, and 33, The Poisson ratios are 0.18, 0.25, and 0.18, shear moduli are 13, 21, and 21. Figure 7 shows that the $y$-displacement error with respect to different discretizations has the same behavior as for isotropic materials (Section 4.3).

4.5. High Aspect-Ratio

To analyze the effect of using elements with high aspect ratio, we repeated the previous experiment for 3 different domains by shrinking the height of the square cross-section from 20, to 4, 2, and 1, while keeping the connectivity identical. This procedure introduces artificial high aspect-ratio elements. To obtain the same anisotropy measure for hexahedral and tetrahedral meshes, we define the aspect ratio between the largest and smallest eigenvalue of the covariance matrix of its vertices. Figure 8 shows the error in the $y$-direction displacement per unit of force for the hexahedral and tetrahedral meshes. The tetrahedral mesh suffers from the low element quality much more than the hexahedral mesh. However, if we regenerate the meshes for the thin domain with the same number of degrees of freedom, but with element quality optimization (last row in Figure 8), the high error of $P_{2}$ disappears and the errors are similar as for $Q_{2}$, as shown in last four rows in Figure 8. For comparison, we also generate a tetrahedral mesh by simply splitting the hexahedra into six tetrahedra. Note that this kind of aspect ratios are extreme and do not appear in any automatically meshed model in our data set (Figure 14).

4.6. 2D Domain with a Hole

Another commonly used test problem is a 2D domain with a hole in the middle. For our experiments we use a square domain of size $200\times 100$ with a hole in the center of radius $20$, the same material model (linear elasticity (4)) and same material parameters $E=210\,000$ and $\nu=0.3$. The experiment consists of applying an opposite in-plane force on the left and right boundary of $100$, that is, stretching the plane horizontally. This problem is obviously ill-posed because of the lack of Dirichlet boundary conditions. We use a standard approach to eliminate the null-space of solutions by exploiting symmetry, and simulating on a quarter of the domain. This leads to a domain with a “corner” cut with two symmetric boundary Dirichlet conditions (displacement is constrained only in the orthogonal direction), a zero Neumann condition, and a Neumann condition corresponding to the original force. We solve this particular benchmark problem on four meshes with different resolutions. Figure 9 shows the von Mises stresses (7) on the top for a triangle mesh and bottom for a quadrilateral mesh, left linear and right quadratic elements. As expected, for a sufficiently dense mesh, all methods converge to similar results. The interesting result is that $Q_{2}$ elements produce visually better results even at really low resolution (first image and second image). In contrast, for linear triangular elements, we need to increase the mesh resolution up to $8\,500$ vertices (last image) for the artifacts to disappear.

This particular problem is also a standard benchmark for incompressible material simulation. We performed the same experiment for a nearly incompressible material: $E=0.1$ and $\nu=0.9999$. Figure 10 shows the norm of the displacement: as for the compressible case, $P_{2}$ and $Q_{2}$ have a similar behavior. Interestingly, for this case, $Q_{1}$ produces a very different solution.

4.7. Nearly Incompressible Material

For the last linear benchmark, we compared the performance of the four discretizations with the material approaching incompressibility. We apply a boundary displacement $[0.2,0]$ on the left and $[-0.2,0]$ on the right of a unit square. We perform a series of experiments in which we keep the Young’s modulus fixed at $0.1$ while changing the Poisson’s ratio from $0.9$ to $0.9999$ ($1$ being the limit of incompressibility in 2D, i.e., area preservation). We compare the standard formulation (4) with a mixed formulation (5) that does not become unstable as $\nu\to 1$. Note that since mixed formulations require different basis degrees for the displacement and the pressure. We performed our experiments using linear pressure bases and quadratic bases for the displacements. We mesh the square with a quad mesh with 4 225 vertices and a tri mesh with 4 229 vertices.

Figure 11 shows the norm of the displacement for this series of experiments. For the nearly incompressible regime (i.e., $\nu=0.9999$) it is remarkable that the quadrilateral element discretization leads to a symmetric and smooth (but incorrect) result for the linear case, while the triangular elements producing an unstable output. The two quadratic discretizations produce visually similar results, close to those obtained with the stable mixed method. The only quantitative difference is that the residual error for the direct solver drops from 1e$-15$ (numerical zero) to 1e$-12$, indicating that the system is close to singular.

4.8. Beam with Torsional Loads

We now compare the solutions for the Neo-Hookean (6) material model for our discretizations. We take a beam with a cross-section $[-10,10]^{2}$ and length 100, $E=200$ $nu=0.35$, fix the bottom part and apply a rotation of $90$ degrees to the top. The rest of the surface is left free. To avoid ambiguities in the rotation we use five steps of incremental loading in the Newton solver. We run the experiment on two sets of tetrahedral and hexahedral meshes. Coarse meshes have $739$ and $750$ while the dense have $50\,000$ and $58\,719$ vertices respectively. The first three images of Figure 12 show that the results are indistinguishable, except for small numerical fluctuations. Similar results are for the plot (Figure 12 last plot) of the rotation angle along a line starting at the point $[9.5,9.5,0]$ parallel to the beam axis. Note that the dense $Q_{2}$ solution required more that 44GB of RAM for the solver, while $P_{2}$ required around 21GB.

The reason for the high memory consumption is the size of the element matrix, which has $(27\times 3)^{=}6\,561$ entries (compared to 144 for $P_{1}$, 576 for $Q_{1}$, and 900 for $P_{2}$). Note that the difference in running time does not come from the number of iterations of the Newton solver: for $P_{1}$, $Q_{1}$, $P_{2}$, $Q_{2}$ we obtained 16, 17, 20, 17 iterations respectively for the coarse mesh, and 18, 16, 17, 18 for the dense mesh.

We have repeated this experiment using quadratic B-spline bases on the coarse mesh. The result is similar to $P_{2}$ and $Q_{2}$, see the inset figure. For this particular example, we measure the solve time of the three discretizations: the spline solve is 3 times faster than $P_{2}$ (0.51s versus 1.50s) and 9 times faster than $Q_{2}$ (0.51s vs 4.62s) while having roughly the same number of iterations: 16. Note that the assembly time (using full integration which could be improved using (Schillinger et al., 2014)) for spline is similar to $Q_{2}$ and is 12 times slower than $P_{2}$ (20.54s versus 1.70s). While using splines on regular grids is natural, the extension to irregular meshes requires the use of T- or U- splines (Beer et al., 2015; Wei et al., 2018), increasing the implementation complexity.

4.9. High Stress

As a final experiment, we run a simulation for an L-shaped domain with the Neo-Hookean material and $E=210000$ $\nu=0.3$. Our goal is to study the differences in the stresses for singular solutions: the concave corner of L will have a stress singularity. We mesh our domains with $14\,155$ vertices for the tetrahedral mesh and $14\,161$ vertices for the hexahedral mesh. We fixed the bottom part of the domain (zero displacement) and rotate the top part by 120 degrees (Dirichlet constraint on the displacement), the rest of the boundary is let free (zero Neumann condition). Figure 13 shows that linear tetrahedral elements underestimate the stress while linear hexahedral elements are somewhat better. The quadratic discretizations are qualitatively similar: the hexahedral mesh does not have the spurious small stress oscillations of $P_{2}$ because the elements are aligned with the mesh, however the price to pay is significant, 17 minutes for $P_{2}$ compared to more than 1.5 hours for $Q_{2}$.

5. Large Dataset

Next, to evaluate the performance of different types of elements on a large diverse set of realistic domains, we compute solutions for the Poisson (1) and linear elasticity (4). We use the method of manufactured solutions (Salari and Knupp, 2000), that is, for an analytically defined solution $u$ we compute the corresponding right-hand side $b$ by plugging it into the PDE. The boundary condition $d$ is obtained by sampling $u$ on the boundary. For the Poisson equation, we use the Franke function (Franke, 1979; Cavoretto et al., 2018)

[TABLE]

while for elasticity

[TABLE]

with Lamè parameters $E=200$ and $\nu=0.35$. In addition to standard tensor product bases for hexahedra, we compare to the popular serendipity bases (Zienkiewicz et al., 2005)[Chapter 6], which have only 20 nodes per element instead of 27.

We use two sources for our data: (1) the Hexalab dataset containing results of 16 state-of-the-art hexahedral meshing techniques (Bracci et al., 2019), (2) the Thingi10k dataset (Zhou and Jacobson, 2016) consisting of triangulated surfaces. For each dataset, we produce a tetrahedral mesh dataset from the surfaces of the hexahedral meshes (generated with MeshGems (Spatial, 2018)) using TetWild (Hu et al., 2018) with a matching number of vertices. Note that since matching the number of vertices is a heuristic process, we discard all meshes where the difference in the number of vertices is larger than 5% of the total number of vertices. To ensure that we are solving a similar problem on the two tessellations we remove meshes whose Hausdorff distance between the surfaces of corresponding hexahedral and tetrahedral meshes differs more than $10^{-3}$ of the diagonal of the bounding box of the hexahedral mesh surfaces. Finally, we discard all meshes whose ratio between boundary and total vertices is more than 30%. Since the Hexalab dataset is small, we opted for doing one step of uniform refinement to increase the number of interior vertices instead of discarding them. In summary, the two datasets are:

1. (1)

237 Hexalab hexahedral meshes and 237 tetrahedral meshes generated with Tetwild. 2. (2)

3 200 hexahedral meshes generated with MeshGems (Spatial, 2018) and 3 200 tetrahedral meshes generated with Tetwild both obtained from the surfaces in the Thingi10k dataset.

For conciseness, we report only the most significant results. Many other metrics (e.g., $H^{1}$-error, the time required to assemble bases, nonzero entries of the matrix, etc.) can be found in the interactive plot.

We remark that, while Tetwild guarantees to produce valid tetrahedral meshes, Meshgems and the methods used in the Hexalab dataset do not provide any guarantee. We observe that out of the 237 Hexalab meshes, 8 (3.4%) contain at least one inverted element (2 from (Livesu et al., 2016) and 6 from (Gao et al., 2016)). For the Thingi10k dataset, Meshgems produces 577 (18.0%) meshes with at least one invalid element. To check if a hexahedron has a negative volume we sample it with $10^{3}$ uniformly spaced samples, evaluate the Jacobian at each point, and mark it as flipped if at least one evaluation is negative. Another important quality measure is the aspect ratio of the elements (Section 4.5). Figure 14 shows that both our datasets contain reasonably well-shaped elements.

All experiments are run on a cluster node with 2 Xeon E5-2690v4 2.6GHz CPUs and 250GB memory, each with max 128GB of reserved memory and 8 threads. For all experiments we use the Hypre (Falgout and Yang, 2002) algebraic multigrid iterative solver and the PolyFEM library for the finite element assembly.

Hexalab

To avoid clutter in the plots we omit the results obtained from meshes with inverted elements leading to plots with 229 points. For the complete statistics see the interactive plot. We first compare the error of the method with respect to the average edge length, Figure 15. We confirm the results of Section 4 for the state-of-the-art hexahedral meshing methods; the accuracy of the solution on a hexahedral and tetrahedral mesh is comparable, in our experiments, for both Poisson and linear elasticity. Figure 16 shows the total and solve time required to reach a certain error, where we draw the same conclusions: the results of the four discretizations are similar. The plots show, as expected, that for a given mesh serendipity elements are faster but less accurate than $Q_{2}$ elements. However this advantage is not consistent enough to change the conclusion related to quadratic tetrahedral elements. Statistics for the individual hexahedral meshing method are available in the interactive plot.

Thingi10k

We repeated the same experiment on $3\,200$ hexahedral meshes generated with MeshGems. For this large dataset it is interesting to note that qualitative behavior of the edge length vs. error curve (Figure 17) is different between hexahedra and tetrahedra: the curve for tetrahedral elements exhibit the expected convergence, while the curve for hexahedra is more flat. This effect comes from the fact that MeshGems is an octree-based method with a tendency to create highly anisotropic meshes. This effect can be mitigated by limiting the difference between the minimal and maximal refinement levels in the octree used to construct the mesh. This leads to more uniform element size, and the results become similar to the results for tetrahedral meshes and the Hexalab dataset, Figure 18.

We also compared running and solve times (Figure 19) and, as expected, serendipity elements are faster than $Q_{2}$ elements but have a larger error. Tetrahedral elements are between the two hexahedral elements: their accuracy is similar to $Q_{2}$ and a running/solving time similar to serendipity.

6. Discussion and conclusions

We presented a large-scale, quantitative study of several common types of finite elements applied to five elliptic PDEs. Our results are consistent on all elliptic PDEs we tried.

We summarize our findings in Figure 20, which allows us to draw the following conclusions for the five elliptic PDEs we considered in our study:

1. (1)

Consistently with well-known observations, $P_{1}$ elements are less efficient (more time spent to obtain a solution with given accuracy) than all other options in all our experiments (Sections 4 and 5). 2. (2)

$Q_{2}$ elements are slightly more accurate than quadratic serendipity elements $\mathrm{SER}_{2}$ but are slightly more expensive for a fixed mesh (Section 5). 3. (3)

$P_{2}$ elements are generally more efficient than $P_{1}$, $Q_{1}$, $Q_{2}$, $\mathrm{SER}_{2}$, that is, we can obtain a given target error in less time, if we can chose the mesh resolution optimal for the desired error level. We were not able to identify any disadvantages for these elements for the range of problems and geometries we have considered (Sections 4 and 5). 4. (4)

Quadratic spline elements $\mathrm{SPLINE}_{2}$ (on a regular lattice) are more efficient than $Q_{2}$ elements (Section 4.8). $\mathrm{SPLINE}_{2}$ are also more efficient compared to $P_{2}$ (3x faster solving time for the same accuracy) but with a much longer assembly time (12 times slower, which could be reduced with more advanced integration techniques (Schillinger et al., 2014)). Their use, however, is restricted by the current meshing technology, as they require meshes with regular grid structure almost everywhere for optimal performance. When these elements are mixed with standard $Q_{2}$ elements to handle general hexahedral meshes with singular vertices and edges  (Schneider et al., 2019a; Wei et al., 2018), their performance advantage is considerably reduced.

For the five elliptic PDEs we considered, unstructured tetrahedral meshes with quadratic Lagrangian basis are a good choice for a “black-box” analysis pipeline: robust tetrahedral meshing algorithms that can process thousands of real-world models exist (Hu et al., 2018), and $p$-refinement can be used to compensate for the rare badly shaped triangles introduced by the meshing algorithms (Schneider et al., 2018).

We leave the extension of this study to non-elliptic PDEs, multiphysics, and collision response as future work. Another important potential extension is the study of bases with orders higher than 2, as is typically the case in IGA setting, or an extension to spectral elements. Another venue for future work is to analyze the impact of the existing different per-element optimizations (e.g., reduced quadrature, hourglass control, special quadrature rules that exploit the tensor-product structure of $Q$ elements, etc.). However, we note that these different optimizations will mostly impact the performance of the assembly and will have little influence on the solve time, which dominates the runtime for sufficiently large problems.

Finally, we acknowledge that the choice of elements is only one of the sources of error in numerical simulations; in realistic scenarios, material models, boundary conditions, or domain shape also play role in the accuracy of a simulation. Extending our study to account for these sources of error is an interesting venue for future work.

Acknowledgements.

We thank the NYU IT High Performance Computing for resources, services, and staff expertise. This work was partially supported by the NSF CAREER award under Grant No. 1652515, the NSF grants OAC-1835712, OIA-1937043, CHS-1908767, CHS-1901091, DMS-1436591, NSERC DGECR-2021-00461 and RGPIN-2021-03707, the SNSF grant P2TIP2_175859, a Sloan Fellowship, a gift from Adobe Research, a gift from nTopology, and a gift from Advanced Micro Devices, Inc.