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

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 packages⁴ 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. Commercial 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. 2005Chapter 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 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. This 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 3,200 hexahedral meshes using the commercial MeshGems-Hexa software Spatial 2018.

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 interactiveplot 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.

Table 1.

View Table

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.

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; and (5) incompressible linear elasticity. We list the corresponding PDEs for completeness.

Let \( \Omega \subset \mathbb {R}^d \), \( d\in \lbrace 2,3\rbrace \) be the domain with boundary \( \partial \Omega \). We aim to solve \( \begin{gather*} \mathcal {F}(x, u, \nabla u, D^2 u) = b,\text{ subject to}\\ u = d \text{ on } \partial \Omega _D\quad \text{and}\quad \nabla u\cdot n = f \text{ on } \partial \Omega _N \end{gather*} \) for the function \( \begin{equation*} u:\Omega \rightarrow \mathbb {R}^n, \end{equation*} \) 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 (1) \( \begin{equation} {\left\lbrace \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}\right.} \end{equation} \)

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 \). (2) \( \begin{equation} {\left\lbrace \begin{array}{ll} -\mu \Delta u + \nabla p = b & \text{on } \Omega \\ -\div {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}\right.} \end{equation} \)

Elasticity. Elasticity PDEs are formulated in terms of the stress tensor \( \sigma [u] \) (which depends on the displacement \( u \)) as (3) \( \begin{equation} {\left\lbrace \begin{array}{ll} -\div {\sigma [u]} = b & \text{on } \Omega \\ u = d & \text{on } \partial \Omega _D \\ \sigma [u] \, n = f & \text{on } \partial \Omega _N. \end{array}\right.} \end{equation} \) 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, (4) \( \begin{equation} \sigma ^L[u] = 2 \mu \epsilon [u]+ \lambda \text{tr} {\epsilon [u]} I \qquad \epsilon [u] = \frac{1}{2} \left(\nabla u^T + \nabla u\right), \end{equation} \) 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: \( \begin{equation*} \mu = \frac{E}{2 (1 + \nu)},\quad \lambda _{\mathrm{3D}} = \frac{E \nu }{(1 + \nu) (1 - 2 \nu)},\quad \mathrm{and}\quad \lambda _{\mathrm{2D}} = \frac{\nu E}{1 - \nu ^2}. \end{equation*} \)

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 \). (5) \( \begin{equation} {\left\lbrace \begin{array}{ll} -\div (2\mu \epsilon [u] + p I) = b & \text{on } \Omega \\ \div 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}\right.} \end{equation} \)

Finally, in the Neo-Hookean material model the stress is a nonlinear function of strain. (6) \( \begin{equation} \sigma [u] = \mu (F[u] - F[u]^{-T}) + \lambda \ln (\det F[u]) F[u]^{-T} \qquad F[u] = \nabla u + I, \end{equation} \) where \( F[u] \) is the deformation gradient.

For elasticity problems, we often use the von Mises stresses (7) \( \begin{equation} \begin{split} 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{split} \end{equation} \) 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.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 three of the four 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.

Fig. 1.

View Figure

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 0 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 four discretizations, using linear and quadratic elements.

Fig. 2.

View Figure

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.

Fig. 3.

View Figure

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.

Fig. 4.

View Figure

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.

Fig. 5.

View Figure

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.

Fig. 6.

View Figure

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).

Fig. 7.

View Figure

4.5 High Aspect-Ratio

To analyze the effect of using elements with high aspect ratio, we repeated the previous experiment for three 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 the last four rows in Figure 8. For comparison, we also generate a tetrahedral mesh by simply splitting the hexahedra into six tetrahedra. Note that these kinds of aspect ratios are extreme and do not appear in any automatically meshed model in our data set (Figure 14).

Fig. 8.

View Figure

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.

Fig. 9.

View Figure

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.

Fig. 10.

View Figure

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 \rightarrow 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.

Fig. 11.

View Figure

4.8 Beam with Torsional Loads

View Figure

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.

Fig. 12.

View Figure

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, and 17 iterations, respectively, for the coarse mesh, and 18, 16, 17, and 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 three times faster than \( P_2 \) (0.51s versus 1.50s) and nine 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=210,\!000 \) \( \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 \).

Fig. 13.

View Figure