Arash Mehraban,Henry Tufo,Stein Sture and Richard Regueiro,★

1Department of Computer Science,University of Colorado Boulder,Boulder,CO,USA

2Department of Civil,Environmental,and Architectural Engineering,University of Colorado Boulder,Boulder,CO,USA

ABSTRACT Higher-order displacement-based finite element methods are useful for simulating bending problems and potentially addressing mesh-locking associated with nearly-incompressible elasticity,yet are computationally expensive.To address the computational expense,the paper presents a matrix-free,displacement-based,higher-order,hexahedral finite element implementation of compressible and nearly-compressible(ν→0.5)linear isotropic elasticity at small strain with p-multigrid preconditioning.The cost,solve time,and scalability of the implementation with respect to strain energy error are investigated for polynomial order p=1,2,3,4 for compressible elasticity,and p=2,3,4 for nearly-incompressible elasticity,on different number of CPU cores for a tube bending problem.In the context of this matrix-free implementation,higher-order polynomials(p=3,4)generally are faster in achieving better accuracy in the solution than lower-order polynomials(p=1,2).However,for a beam bending simulation with stress concentration (singularity),it is demonstrated that higher-order finite elements do not improve the spatial order of convergence,even though accuracy is improved.

KEYWORDS Matrix-free;higher-order;finite element;parallel;linear elasticity;multigrid solvers;unstructured meshes

1 Introduction

Modeling bending deformation and associated stress of compressible and nearly-incompressible linear isotropic elastic materials via the standard linear interpolation displacement-based Finite Element Method (FEM) may lead to inaccurate stress calculations (depending on element sizeh) as well as potential mesh-locking [1,2].With respect to accuracy,eitherh-orp-refinement are adopted,but oftentimes not the two together because of the associated high computational cost.With respect to mesh-locking,several FEM approaches involving mixed formulations have been proposed to overcome volumetric strain locking of linear finite elements [3–6].Most such mixed methods for linear elasticity at small strain solve for two fields (displacement and pressure),which are discretized separately via an additive split of the deviatoric and volumetric parts of the small strain tensor [7,8].As an alternative to,and oftentimes in conjunction with,mixed formulations,higher-order finite elements have been adopted to attempt to overcome mesh locking [9].However,assembling a global Jacobian matrix based on higher-order finite element interpolation functions leads to expensive solve times as a result of higher memory requirements [10].Matrixfree formulations [11,12] provide the benefits of higher-order finite element methods but with more efficient memory usage.Concurrently,using tensor-product-basis evaluation further improves the performance of higher-order methods by supplying favorable storage and work estimates ofO(pd) andO(pd+1) perp-order element,respectively,for discretizations in Rd[11],whereddenotes spatial dimension.Furthermore,tensor-product-based operator evaluations can be cast as matrixmatrix products.Iterative solvers,such as Krylov subspace methods,only require the result of matrix-vector products rather than an assembled matrix,which allows these methods to take advantage of the computational efficiencies offered by matrix-free formulations.Therefore,storage of large,sparse matrices is avoided with iterative solvers.

Krylov subspace methods require preconditioning [13] to efficiently and stably solve largescale elliptic problems such as 3D linear isotropic elasticity within the static balance of linear momentum equations.Geometric multigrid is a robust preconditioning scheme that is an appealing choice for structured meshes,and has been considered in numerous studies using finite difference and finite element methods with tensor-product-bases implemented on CPUs and GPUs[14–16].p-multigrid,developed by Rønquist et al.[17],is a version of geometric multigrid based on coarsening by decreasing the basis order in higher-order or spectral finite elements rather than coarsening by aggregating elements.As a result,p-multigrid is a natural approach for problems on unstructured meshes,whereby spatial convergence only weakly depends on polynomial orderp,if properly implemented.Computational cost of the discretization scheme may be affected by polynomial orderpand element sizeh,among other factors [18,19],such as adaptiver-refinement of the mesh [20].

In terms of engineering applications,compressible linear isotropic elasticity (e.g.,ν=0.3) is appropriate for modeling metals in their elastic regime,such as in serviceability studies to analyze the stress intensity factor at a crack tip to determine if the crack will propagate and thus estimate the life cycle of a metallic component.Nearly-incompressible linear isotropic elasticity (ν→0.5)is appropriate for modeling small deformations (small strains and small rotations) of rubber-like materials such as solid rock propellant binders [21].When metals or rubber-like materials are loaded beyond their small strain linear elastic limit,then large deformation hyperelasticity or hyper-elasto-plasticity constitutive models are needed,which are beyond the scope of this paper.

In this work,we apply a previously-developed matrix-free,higher-order,FEM for solving three-dimensional (3D) linear isotropic elasticity withp-multigrid preconditioning [22,23] to (i) a tube geometry subjected to method of manufactured solutions (MMS),(ii) a tube bending problem for compressible (ν=0.3) and nearly-incompressible (ν=0.499999) elasticity and polynomial orderp=1,2,3,4,and (iii) a beam bending problem under body force loading for which the beam is clamped on both ends.Simulations (ii) and (iii)—but simulation (iii) in particular—cause a stress concentration (singularity) in the clamped regions,where the higher-order FEM may not converge optimally as expected [22].Therefore,we investigate order of convergence of the matrixfree higher-order FEM for compressible elasticity when there is a stress concentration.We aim to answer the following questions:(1)With the matrix-free parallel implementation,what combinationof mesh-refinement h and polynomial order p leads to a certain strain energy error that requires the least computational cost for compressible and nearly-incompressible elasticity using the standard displacement-based FEM?(2)What combination of mesh-refinement h and polynomial order p achieves a certain strain energy error for least amount of simulation time?(3)Does the higher-order,displacementbased FEM remain competitive as Poisson’s ratio approaches the nearly-incompressible limit(ν→0.5)?(4)Does the higher-order FEM converge as expected when there is a stress concentration present?The parallel computational toolkit PETSc [24] is employed for the linear solver as well as Algebraic MultiGrid (AMG) preconditioner for the coarse-grid solve,while libCEED [25] is used to perform efficient tensor-product-basis evaluation.

The rest of the paper is organized as follows.In Section 2,the constitutive model for compressible 3D linear isotropic elasticity and the variational form of static balance of linear momentum are presented.In Section 3,a preconditioning scheme used to accelerate convergence is briefly discussed.In Section 4,numerical results are presented,and in Section 5,observations and conclusions are provided.

2 Compressible Linear Isotropic Elasticity

For linear isotropic elastic materials,we consider the stored strain energy function as,

Therefore,its stress-strain relationship is given by its derivative with respect to strain as,

whereσand∈are stress and strain tensors,respectively;andλandμare the Lamé parameters.The small strain tensor is,

The strong form of the static balance of linear momentum is given by the following,

2.1 Residual Evaluation

Discretization of the variational Eq.(5) can be written to facilitate matrix-free evaluation [10].The residual in discrete form can be expressed as,

where Neand Beiare evaluations of the finite element shape functions and their derivatives at the quadrature points inx,y,andzdirections,Eeis the elementerestriction operator that separates Degrees of Freedom (DoF) based on the elements to which they belong,and represents pointwise function evaluation.f0andf1are determined from the constitutive model and its tangent,whereue=Ne(Eed)andanddis the total nodal displacement vector.The basis operators are represented as Kronecker products,

2.2 Jacobian Evaluation

Similar to a residual evaluation,the action of the Jacobian can be computed using the notation proposed by [10,26],as,

where

In the small strain case,for Eq.(2),f0is not a function ofuor ∇u.Therefore,its derivative with respect touor ∇uis zero,(i.e.,f0,0=0 andf0,1=0).On the other hand,f1is a function of ∇u,but it is not a function ofu.Therefore,f1,0=0 andf1,1=∂σ/∂∈;due to linearity,∂σ/∂∈=λI⊗I+2μI,or equivalently,f1,1(∇δu)=(∂σ/∂∈):∇δu=λtr[∈(δu)]+2μ∈(δu),which is equivalent to (2) applied to an incrementδuwithin a nonlinear iterative solver.While these minor simplifications are possible for linear problems in the present work,our implementation solves the problem as though it were nonlinear,and we will continue using the corresponding terminology.

3 p-Multigrid Preconditioning of Linear Isotropic Elasticity

Using the formulation in Eq.(8),we can compute the action of global Jacobian matrix on with arbitrary user defined polynomial orderp.An iterative linear equation solver is required with matrix-free operators,which necessitates preconditioning,especially at higher-order.With unstructured meshes,a natural hierarchy of grids does not exist,soh-multigrid can be difficult to implement.Algebraic MultiGrid (AMG) is suitable for low order meshes for which the Jacobian matrix can be assembled.However,assembly of this matrix is prohibitively expensive for higherorder element meshes [27].Therefore,we use geometric multigrid withpcoarsening while adopting AMG as the coarse grid solver.A Chebyshev polynomial smoother based upon the operator diagonal [28] is called in the multigrid cycle.

Inp-multigrid,grid transfer operations increase or decrease the polynomial order of the element basis functions,and these operations can be implemented in a matrix-free fashion viaf0,0of (8) with suitable basis evaluatorsN.The coarse-to-fine (ctof) basis operation,Nctof,interpolates the DoFs on the nodes of a coarse grid element to the nodes of a fine grid element (N27×8forQ1prolongation toQ2,for example,whereQ1denotesp=1 andQ2,p=2).The corresponding coarse and fine grid element restriction operators,Ee,candEe,f,are used in the grid transfer operators.The operatorcorrectly computes the interior degrees of freedom but over-counts nodes on the element facets.We can count the multiplicity of each node on the fine grid by applying the transpose fine grid restriction to the unit vector,Thus,thep-multigrid prolongation operator is given by,

and then thep-multigrid restriction operator is given byR=PT.

4 Numerical Examples

Using the 3D linear isotropic elasticity model and matrix-free,higher-order FEM presented in Sections 2 and 3,we present numerical results based on (i) method of manufactured solutions(MMS) applied to a tube geometry,(ii) a tube bending simulation for compressible and nearlyincompressible elasticity on unstructured meshes,and (iii) a beam bending simulation with stress concentration on structured meshes.Polynomial ordersp=1,2,3,4 are applied for a range of meshes in the compressible case,and polynomial ordersp=2,3,4 are applied for the same range of meshes in the nearly-incompressible case.For the compressible case,Poisson’s ratioν=0.3 and Young’s modulusE=69×109Pa,that correspond to aluminum.For the nearly-incompressible case,ν=0.499999 andE=0.1×109Pa,that correspond to rubber-like materials in the small elastic strain regime.

The parallel platform used to run all FE simulations with higher-order polynomial element meshes is an Intel Xeon E5-2680 v3 @2.50 GHz (2 CPUs/node,24 cores/node) and 113 GB RAM per compute node [29].The PETSc toolkit v3.14 is called for mesh management,domain decomposition,parallel assembly operations,and communication over all MPI processes.The matrix-free Jacobian operator is implemented using libCEED v0.8 and PETSc for Eq.(8).Preconditioning withp-multigrid and AMG is enabled through PETSc’s multigrid interface with grid-transfer and operator application at each level implemented in libCEED.Vectorized tensor-product operations for 8-element batches in the matrix-free operator evaluation on AVX-2 extensions of×86 instruction sets for each local processor are conducted in libCEED.For simulations conducted with MMS and tube bending,an unstructured mesh consisting of 400 Hex8 elements (trilinear) is generated in the Trelis [30] software.The mesh is then refined 5 times to generate 3,200,10,800,25,600,50,000,and 86,400 element meshes,respectively.

4.1 Achieving L2 Strain Energy Error via h-and p-Refinement in Tube Meshes Subjected to MMS

Considering the tube geometry in Fig.1 with method of manufactured solutions (MMS)-generated displacement boundary conditions,the first four tube unstructured mesh refinements with 3,200,10,800,25,600,and 50,000 elements are used in the FE simulations for the compressible and nearly-incompressible cases.For the MMS,we produce a contrived,but smooth,right hand sideρbased on=[u1,u2,u3]Twithu1=e2xsin(3y) cos(4z)/108,u2=e3xsin(4y) cos(2z)/108,andu3=e4xsin(2y) cos(3z)/108.Every problem is run once to determine theL2global strain energyΦherror from the MMS.Tables 1a and 1b summarize the size of each problem in terms of degrees of freedom (“#DoF”) based onh-andp-refinement and the correspondingL2strain energy error based on MMS for the compressible case.Tables 2a and 2b summarize the size of each problem in terms of DoF based onh-andp-refinement and the correspondingL2strain energy error based on MMS for the nearly-incompressible case.Considering the smoothness of the MMS,and relatively high number of DoF for even the coarsest mesh (3009 elements),the strain energy error is already quite small prior to refinement.Although,forp-refinement,we see for the compressible case some convergence,while for the nearly-incompressible case,the error is small and stays small (but not as small as the compressible case,10-10vs.10-7) indicating stability for smooth deformations via MMS.

Table 1:L2 error (using MMS on tube geometry) for compressible elasticity (ν=0.3) with polynomial orders,p=1,2,3,4

Table 2:L2 error (using MMS on tube geometry) for nearly-incompressible elasticity(ν=0.499999) with polynomial orders,p=1,2,3,4

Figure 1:Deformed hexahedral mesh for the cylindrical tube bending problem.The length of the tube is 100 mm,with circular cross-section with inner diameter 15 mm and outer diameter 20 mm.The left wall is fixed in x,y,and z displacements,while the right wall is displaced in the negative y direction by 3 mm

4.2 Achieving Reduced L2 Strain Energy Error via h-and p-Refinement in Tube Meshes Subjected to Bending

Tube meshes subjected to bending are simulated to study cost and spatial convergence of the matrix-free,higher-order FEM forh-refinement andp-refinement based on strain energyΦhfor compressible and nearly-incompressible elasticity.The length of the tube is 100 mm,with circular cross-section of inner diameter 15 mm and outer diameter 20 mm.A coarse mesh representation of the tube is shown in Fig.1.The left side of the tube is clamped by placing a zero displacement boundary condition on the surface,while a second displacement boundary condition is imposed on the right surface of the tube to bend it 3 mm in the negativeydirection.Fig.1 shows the deformed mesh after bending.All tube meshes are simulated for compressible and nearly-incompressible elasticity.Tables 3 and 4 summarize the problem sizes in terms of degrees of freedom (“#DoF”) for each refinement in conjunction with different polynomial orderspfor compressible and nearly-incompressible elasticity,respectively.The size of problem (“#DoF”)increases when a finer meshhor a higher-order polynomialpis used.As the FE mesh is refined inh,the norm of the strain energy computed from the FE solution converges to the exact norm of strain energy [31].Therefore,we use the discrete elastic strain energyΦhto compute relative error in the FE simulations with respect to smaller size problems.Achieving a certain strain energy error in the FE solution is problem-dependent.For example,problems with smooth and nonsmooth solutions behave differently with mesh refinementhwhen used with different polynomial ordersp[18].Therefore,with mesh refinement through decreasingh,we treat the strain energy from the finest mesh as the “exact” solution in the computation of relativeL2error for each polynomial orderpto determine the effect ofh-refinement in achieving a certain error.However,withp-refinement,we treat the strain energy from the highest polynomial order (p=4) to compute the relativeL2strain energy error for eachp-refinement.Tables 3a and 3b summarize the relativeL2error based on mesh refinementhand using higher-order polynomialspfor the compressible tube bending FE simulations,respectively.In addition,depending on the problem,each simulation is run with 1,6,12,24,48,96,192,384,and 768 CPU cores for the compressible case.Each simulation is run 3 times to reduce noise in the performance timing for a total of 297 simulations for the compressible case.

Table 3:L2 error based on strain energy Φh for compressible (ν=0.3) elastic tube bending

Table 4:L2 error based on strain energy Φh for nearly-incompressible (ν=0.499999) elastic tube bending

(Continued)

Table 3(continued)p #Refine #DoF Strain Energy L2 Error(a) L2 error based on h-refinement 4 1 686880 3.473070e+00 7.7888e-03 4 2 2237040 3.490239e+00 2.8839e-03 4 3 5206080 3.496170e+00 1.1896e-03 4 4 10054800 3.498879e+00 4.1542e-04 4 5 17244000 3.500334e+00 0.0000e+00#Refine p #DoF Strain energy L2 Error(b) L2 error based on p-refinement 1 1 14040 3.600109e+00 3.6578e-02 1 2 94800 3.483284e+00 2.9409e-03 1 3 299880 3.475215e+00 6.1749e-04 1 4 686880 3.473070e+00 0.0000e+00 2 1 42480 3.554248e+00 1.8339e-02 2 2 299880 3.495535e+00 1.5174e-03 2 3 966600 3.491421e+00 3.3857e-04 2 4 2237040 3.490239e+00 0.0000e+00 3 1 94800 3.535690e+00 1.1304e-02 3 2 686880 3.499542e+00 9.6448e-04 3 3 2237040 3.496988e+00 2.3394e-04 3 4 5206080 3.496170e+00 0.0000e+00 4 1 178200 3.526167e+00 7.7990e-03 4 2 1313400 3.501305e+00 6.9319e-04 4 3 4305600 3.499495e+00 1.7597e-04 4 4 10054800 3.498879e+00 0.0000e+00 5 1 299880 3.520558e+00 5.7779e-03 5 2 2237040 3.502210e+00 5.3612e-04 5 3 7366680 3.500819e+00 1.3882e-04 5 4 17244000 3.500334e+00 0.0000e+00

To address the first question raised in the Introduction with regard to minimal computational cost for a certain strain energy error in the FE solution for compressible elasticity,we compute the minimum amount of CPU work to reach a certain error based on mesh refinements and higher-order polynomials.To do this,for each polynomial orderpwe accumulate the time spent to reach the error.Then we determine the minimum time in the list and multiply it by the number of processorsnpthat are executed for the minimum time.We generatePareto optimaldiagrams to present the results of the study for compressible and nearly-incompressible elasticity separately,where a point on the diagram is calledPareto optimalif strain energy error cannot be decreased(moving down the ordinate (yaxis)) without increasing time (or cost) by moving to the right on the abscissa.The set of Pareto optimal configurations is known as thePareto front.

Figs.2a and 2b present theL2strain energy error based on mesh refinementvs.cost,and theL2error based on higher-order polynomialsvs.cost,for compressible elasticity.Figs.2a and 2b show that the minimal cost to achieve a certain error occurs in the lower left region of the diagram,indicating that higher orderpis the most cost-effective for the compressible case.Each horizontal series of the same color represents a strong scaling study at fixed mesh refinementhand polynomial orderp,but with changing number of processorsnp,with perfect strong scaling occurring when all color dots are collocated.The more expensive simulations tend to exhibit better strong scaling because they have more work over which to amortize the inherent communication costs,while smaller models (fewer DoFs) are more cost-efficient to run on a single core.The low orderp=1 cases are increasingly far from the Pareto front.In addition,solely calling larger number of CPU cores is less cost-efficient,although are faster with respect to solve time up to the point where the parallel overhead communication begin to dominate.

Figure 2:Error vs.cost for compressible elastic tube bending.The Pareto optimal configurations occur in the lower left region of the plot indicating a higher-order polynomial p is more cost-efficient for a certain error.(a) Error based on h-refinement vs.cost.(b) Error based on p-refinement vs.cost

To answer the second question raised in the Introduction regarding the least amount of time to achieve a certain strain energy error given enough computational resources,we provide a second set of Pareto optimal diagrams in Figs.3a and 3b for the compressible elasticity tube bending problem.The Pareto optimal configurations are toward the lower left region of the plot,whereby uponp-refinement,p=3 exhibits faster solve times for certain error,while uponh-refinement,p=2 delivers the fastest solve times for certain error,a result that was unexpected and requires further investigation.Each horizontal series of the same color represents strong scaling of a given mesh refinementhand polynomial orderp.Larger number of CPU cores when enough work is amortized tend to reach a certain error in the solution faster with higher order polynomials.Similar to Pareto optimal for the cost,strong scaling occurs when the same color larger dots overlay smaller dots.

Figure 3:Error vs.time for compressible elastic tube bending.The Pareto optimal configurations occur toward the lower left region of the plot.(a) Error based on h-refinement vs.time.(b) Error based on p-refinement vs.time

We notice that depending on the error,increasing the polynomial orderpwhile fixing mesh size (hconstant) may provide the fastest time to solution.For example,for strain energy error between 10-3and 10-4,polynomial orderp=3 may be faster than lower order polynomials,but for error greater than 10-3,polynomial orderp=2 may be faster.This could be due to the number of DoFs for different polynomial orders when run with the same size mesh.In addition,we expected that higher order polynomialsp=3,4 would be more cost-efficient than lower order polynomials.However,we observed that polynomial orderp=2 is more cost-efficient than higher order polynomialsp=3,4 for this compressible tube bending problem.

For nearly-incompressible elasticity,we use the same mesheshand polynomial orderp.The problem sizes in terms of DoFs andL2strain energy error based on mesh refinement and polynomial orders are represented in Tables 4a and 4b.However,we instead use 24,48,96,192,384,768,and 1536 CPU cores since nearly-incompressible elasticity is more computationally intensive for matrix-free solution.This resulted in 228 FE simulations.Figs.4a and 4b represent Pareto optimal diagrams for theL2error based on mesh refinementvs.cost,and theL2error based on higher-order polynomialsvs.cost,respectively.Pareto optimal consists of polynomial ordersp=2,3 whenL2error is computed based onh-refinement,while polynomial orderp=3 is more cost-effective in achieving more accuracy in the solution.Similarly,the Pareto optimal diagrams are provided forL2error based on mesh refinementh vs.solve time,andL2error based on polynomial orderp vs.solve time in Figs.5a and 5b,respectively.The Pareto front consists ofp=2,3,4 in that order.In addition,larger number of CPU cores provide faster solutions for a certain strain energy error.

Figure 4:Error vs.cost for nearly-incompressible elastic tube bending.The Pareto optimal configurations occur toward the left region of the plot.Perfect strong scaling occurs when same color dots are collocated.(a) Error based on h-refinement vs.cost.(b) Error based on p-refinement vs.cost

Figure 5:Error vs.time for nearly-incompressible elastic tube bending.The Pareto optimal configurations occur toward the left region of the plot.Perfect strong scaling occurs when same color dots are collocated.(a) Error based on h-refinement vs.time.(b) Error based on p-refinement vs.time

Next,we examine the relative parallel scalability of the matrix-free FE implementation in terms of throughputvs.solve time forp=1,2,3,4 in the compressible case,andp=2,3,4 in the nearly-incompressible case.Using the tube geometry,new meshes are generated in Trelis.A mesh with 1,638,400 elements is used with polynomial orderp=1.Meshes with 204,800,60,500,and 25,600 elements are used with polynomial orderp=2,3,4,respectively.These mesh sizes in conjunction with polynomial ordersp=1,2,3,4 produce approximately 5.2M DoFs.This problem size is conveniently analyzed in the Pareto diagrams,which roughly occupies one compute node (24 CPU cores) on the parallel platform.For the compressible case,depending on the number of elements,12,24,48,96,192,384,and 768 CPU cores are executed where every problem is run three times for a total of 75 FE simulations.For the nearly-incompressible case,192,384,768,and 1536 CPU cores are called.This resulted in 42 FE simulations for which each problem is run three times.Throughput is computed as millions of DoFs per second per CPU core for each FE simulation.The result of throughputvs.solve time for the compressible and nearly-incompressible cases are presented in Figs.6a and 6b,respectively.We notice for the compressible case,the implementation scales up to 200–250 CPU cores (Fig.7a) with polynomial orders 3 and 4.The sub-optimal performance in the code pertains to lack of work amortized to each CPU core for the problem size chosen.However,for the nearly-incompressible case,which is more computationally-intensive compared to the compressible case,the implementation scales up to 800–1000 CPU cores (Fig.7b) for polynomial orders 2–4 for the problems sizes selected.

Figure 6:Throughput vs.solve time for compressible and nearly-incompressible elasticity problems.For the problem size with 5.2 M DoFs across polynomial orders p=1,2,3,4.(a) Throughput vs.solve time.(b) Throughput vs.solve time

4.3 Beam Bending under Body Force

In this section,we consider a solid square-cross-section aluminum beam subject to a body force in the negativeydirection simulated withp=1,2,3,4.Fig.8 shows a coarse structured mesh of the beam.The length of the beam is 5 m,and each side of the square beam crosssection is 0.25 m.The beam is clamped on both ends.A body force of 200N/m3is applied in the negativeydirection on the beam,which causes the beam to bend 9.1×10-4mm in the negativeydirection at the middle of the beam.Fig.8 shows the deformed mesh of the beam.Similar to the tube bending simulation,the clamped BCs will generate stress concentration(singularity) in the clamped regions.When a stress concentration exists,it can be shown that the order of convergence of the FE solution becomes independent of polynomial orderp[32].That is,the order of convergence of the FE solution cannot be improved by using higher-order finite elements.We investigate this phenomenon by conductingh-andp-refinement with our matrix-free FE implementation in libCEED/PETSc.

Figure 7:Throughput vs.number of processors for compressible and nearly-incompressible elasticity problems.For the problem size with 5.2 M DoFs across polynomial orders p=1,2,3,4,the implementation remains scalable up to 200–250 CPU cores in the compressible case with polynomial order 3 and 4,while it is scalable up to 800–1000 CPU cores in the more computationallyintensive nearly-incompressible case.(a) Throughput vs.solve time.(b) Throughput vs.solve time

Figure 8:Deformed hexahedral beam structured mesh 5 m in length and square cross-section with side 0.25 m.Beam bent by 200 N/m3 body force in the negative y direction,while the left and right ends are clamped.The body force causes a 9.1×10-4 mm displacement in the middle of the beam

Five different refinements of the Hex8 beam mesh are generated with the Trelis software for 558,4,464,15,066,35,712,and 69,750 elements.The finest mesh with 69,750 elements acts as the “exact” solution with respect to strain energy calculation.A total of 60 simulations are conducted where each simulation is repeated 3 times.AnL2strain energy error based onh-refinement per polynomial order is computed forp=1,2,3,4.In addition,anL2strain energy error based on higher-order polynomialpfor each mesh refinementhis computed.Tables 5a and 5b summarize the refinement (“#Refine”),polynomial degreep(“deg”),number of degrees of freedom (#DoF),strain energy,and theL2error.The element size for each mesh is computed to beh=0.1428,0.0714,0.0476,0.0357 m for smallest to largest number of elements mesh,respectively.Fig.9 represents theh-refinement order of convergence for polynomial orders 1 through 4.The calculated slopes (order of convergence) of the lines that are best fit in the least square sense for polynomial ordersp=1,2,3,4 withh-refinement are 2.43,2.48,2.38,and 2.31,respectively.Therefore,we conclude that for this problem the order of convergence of the FE solution on structured meshes with stress concentration (singularity) cannot be improved by using higher-order polynomials.Thus,the lack of optimal spatial convergence with regard to higherorder FE when there is a stress concentration is confirmed,as opposed to when there is no stress concentration (or singularity) [22].

Table 5:L2 error based on strain energy Φh for compressible elasticity (ν=0.3) beam bending under body force for polynomial orders,p=1,2,3,4

(Continued)

Table 5(continued)#Refine p #DoF Strain Energy L2 Error(b) L2 error based on p-refinement 0 1 2928 1.585403e-05 3.1509e-02 0 2 18081 1.633363e-05 2.2117e-03 0 3 55500 1.636293e-05 4.2167e-04 0 4 125229 1.636983e-05 0.0000e+00 1 1 18081 1.621649e-05 9.5842e-03 1 2 125229 1.636326e-05 6.2028e-04 1 3 401793 1.637133e-05 1.2775e-04 1 4 928125 1.637342e-05 0.0000e+00 2 1 55500 1.629746e-05 4.6936e-03 2 2 401793 1.636923e-05 3.1068e-04 2 3 1310064 1.637321e-05 6.7510e-05 2 4 3051501 1.637431e-05 0.0000e+00 3 1 125229 1.632844e-05 2.8228e-03 3 2 928125 1.637149e-05 1.9338e-04 3 3 3051501 1.637394e-05 4.3745e-05 3 4 7138173 1.637466e-05 0.0000e+00 4 1 237312 1.634363e-05 1.9060e-03 4 2 1784577 1.637263e-05 1.3509e-04 4 3 5897292 1.637432e-05 3.1825e-05 4 4 13830957 1.637484e-05 0.0000e+00

Figure 9:h-refinement order of convergence plot for clamped beam bending under body force with polynomials p=1,2,3,4

5 Conclusions

The paper investigated the use of displacement-based higher-order finite elements for compressible and nearly-incompressible linear isotropic elasticity using a matrix-free implementation withp-multigrid preconditioning.It was observed that varying the polynomial orderpand mesh sizehaffects the cost and time it takes to achieve a certain strain energy error for Poisson’s ratiosν=0.3,0.499999 in a tube bending problem.Increasing the polynomial orderpwhile fixing the mesh size (hconstant) may provide the fastest time to converge to a certain strain energy error between 10-3and 10-4.However,for strain energy error greater than 10-3,lowerorder polynomials may achieve the solution faster.When there is a stress concentration present,higher-order polynomials do not improve the order of convergence in the finite element solution.In addition,polynomial orderp=2 proved to be more efficient thanp=3,4 for the compressible elastic tube bending problem,which also has a stress concentration (singularity in the stress solution).Each mesh refinement was run with polynomial ordersp=1,2,3,4 for compressible elasticity,and polynomial ordersp=2,3,4 for nearly-incompressible elasticity.Using the same mesh refinementshwith higher-order polynomials (p=3,4) produce a much larger problem size (more DoFs) than when the same meshes are used with polynomial order 2.Therefore,polynomial ordersp=3,4 appear to be more costly than polynomial orderp=2.An exhaustive search using different polynomial orderspand mesh sizeshmay determine if polynomial orderp=2 is truly more cost efficient than polynomial ordersp=3,4,which will be considered as future work.In addition,for nearly-incompressible elasticity (ν=0.499999),the tube bending problem could not achieve an error better than 10-1based on discrete global strain energyΦhcomputation because of the stress concentration;whereas a smooth method of manufactured solution (MMS) tube problem achieved strain energy error on the order of 10-7,illustrating the difficulty in convergence forh-andp-refinement for near-incompressibility when a stress concentration is present.Therefore,implementing a mixed formulation for handlingν→0.5 appears to be expedient with respect to spatial convergence than solely a displacementbased,higher-order finite element method.Scalability of the matrix-free FEM for compressible elasticity is feasible up to a few hundred CPU cores,while scalability improves for the more computationally-expensive nearly-incompressible case to one thousand CPU cores for the problem sizes chosen in this study.Higher-order finite elements generally still out perform lower-order finite elements within the matrix-free implementation presented in the paper,for which the cost per degree of freedom decreases with increasingp,due to more structured computation and more efficient quadrature.

Notation:Boldface denotes vectors and tensors in symbolic form.Unless otherwise indicated,all vector and tensor products in symbolic form are assumed to be inner products,such asνν=vivi,(ab)ik=aijbjkand(a):(b)=aijbij,where repeated indices denote a sum over those indices.Cartesian coordinates are assumed.The symbol tr(·) is the trace operator,such that tr(σ)=σii.The symbol I is the unit tensor,i.e.,Iij=δij,whereδijis the Kronecker delta operator,and I is the fourth order identity tensor.prefers to polynomial basis order,hto characteristic element size,dto spatial dimension,andnpto number of processors.

Funding Statement:The research relied on computational resources [29] provided by the University of Colorado Boulder Research Computing Group,which is supported by the National Science Foundation (Awards ACI-1532235 and ACI-1532236),University of Colorado Boulder,and Colorado State University.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.