Applied Mathematics and Computation A posteriori discontinuous Galerkin error estimator for linear elasticity

This paper presents for the ﬁrst time the derivation of an hp a posteriori error estimator for the symmetric interior penalty discontinuous Galerkin ﬁnite element method for linear elastic analysis. Any combination of Neumann and Dirichlet boundary conditions are ad- missible in the formulation, including applying Neumann and Dirichlet on different components on the same region of the boundary. Therefore, the error estimator is applicable to a variety of physical problems. The error estimator is incorporated into an hp -adaptive ﬁnite element solver and veriﬁed against smooth and non-smooth problems with closed- form analytical solutions, as well as, being demonstrated on a non-smooth problem with complex boundary conditions. The hp -adaptive ﬁnite element analyses achieve exponential rates of convergence. The performances of the hp -adaptive scheme are contrasted against uniform and adaptive h reﬁnement. This paper provides a complete framework for adap- tivity in the symmetric interior penalty discontinuous Galerkin ﬁnite element method for linear elastic analysis.


Introduction
The discretization of partial differential equations for numerical computation facilitates the solution of physical systems, however, it introduces approximation errors. Error estimators are essential for the numerical analysis of boundary value problems as they allow for the assessment of the accuracy of a simulation in the absence of analytical solutions [1] .
In this work, we present a new error estimator for linear elasticity based on a Discontinuous Galerkin (DG) finite element method. DG methods were introduced in the early 1970s as a way to numerically solve first-order hyperbolic problems [2] . More recently, these methods have been applied to elliptic problems [3][4][5] . A variety of DG methods have been developed in the decades [6] , among them is the subclass of interior penalty DG methods which are stable and, in the authors' opinion, easy to implement. For these reasons the method used in this work is the symmetric version of the interior penalty DG method.
Error estimators for elliptic problems have been an important topic in the last few years. In [7] the first residual-based energy-norm error estimator for the symmetric interior penalty discontinuous Galerkin finite element method is presented. In the following years, inspired by this pioneering work, error estimators have been derived to solve a variety of problems (see for example [8][9][10][11] ). The present work is meant to fill the need of an hp residual-based energy-norm error estimator for DG methods for linear elasticity.
The layout of the paper is as follows. Section 2 introduces the linear elastic model problem including pure Neumann/Dirichlet and mixed Neumann/Dirichlet boundary conditions and Section 3 presents the weak formulation for this problem. The symmetric interior penalty discontinuous Galerkin finite element method is introduced in Section 4 , and Section 5 provides the a priori convergence results of the method. Sections 6 and 7 present the reliability and the efficiency of the error estimator. The analysis in Sections 6 and 7 is only presented for the two dimensional case, but it holds also for the three dimensional case. The restriction to the two dimensional case is to keep the paper more readable. For the same reason we only consider triangular elements, but the analysis is also applicable to affine quadrilateral elements and to tetrahedral and affine hexahedral elements in three dimensions. Numerical examples verifying the error estimator and its implementation within an hp -adaptive solver are presented in Section 8 before conclusions are drawn in Section 9 .

Model problem
The model problem considered in this paper is linear elasticity with several kinds of boundary conditions: let be a bounded polygonal domain in R 2 with ∂ = D ∪ N ∪ T , where D , N and T are disjoint sets, and let u the solution  (1) imposes that the tangential component of the traction is zero.
In order to ensure that the model problem (1) has a unique solution, the rigid motions have to be excluded. In order to do that lets introduce the function space S := [ H 1 ( )] 2 \ R, where R is the space containing all rigid motions, clearly rigid motions are not in S. In view of this we assume that u ∈ S. The lack of uniqueness of the solution can be traced to the fact that the kernel of the strain operator ( · ) contains only rigid motions [12] .

Remark 1.
One or more sets among D , N and T can be empty. We assume that all the considered combinations of boundary conditions ensure the uniqueness of the solution of problem (1) up to rigid motions. A way to filter out the rigid motions in simulations is presented in [13] . This technique is used in Section 8.3 to solve the crack in a plate problem.
We define the strain tensor for a displacement v as (v ) i j := 1 2 (∇ j v i + ∇ i v j ) and we restrict the choice of material properties such that we have that σ(v ) = D (v ) where the matrix D ∈ R 3 ×3 ×3 ×3 is symmetric and invertible implying that there are two positive constants D min and D max such that It is straightforward to see that on S there are two positive constants c and C such that with · 0, the L 2 norm. Also, thanks to the fact that D is invertible, there are two positive constants c σ and C σ such that (4) Theorem 2.1 (Coercivity in the continuous case) . There is a positive constant c CG such that for any function v ∈ S we have: Proof. The statement comes easily from (3) and (4) .

Continuous weak formulation
In this section we introduce the continuous weak formulation of problem (1) , this formulation is only used in the analysis for the error estimator and it is not implemented in the code. Due to the variety of boundary conditions in (1) , we introduce different spaces for the trial and test functions: with n the unit vector perpendicular to the boundary of and pointing outward. Thus, the weak form of problem (1) reads as follows: find u ∈ S CG 1 such that where the bilinear form a CG ( ·, ·) and the linear form l CG ( · ) are The definitions in (7) can be derived from the strong problem (1) applying integration by parts and then applying the boundary conditions.

Symmetric interior penalty discontinuous Galerkin method
In this section, we introduce our DG method to solve problem (1) on . Throughout, we assume that the computational domain can be partitioned into a shape regular mesh T of triangular elements and we denote with K a generic element of T . Also we assume that the elements are affine, i.e. for each element K there exists an affine map between the reference element ˆ K and the physical element K . We allow for a maximum of one hanging node per edge and we denote E (T ) and E int (T ) ⊂ E (T ) the set of all edges of the mesh T and the subset of all interior edges, respectively, and by E BC (T ) ⊂ E (T ) the subset of all boundary edges. Furthermore, the set E BC (T ) is partitioned in the three subsets E D (T ) , E N (T ) and E T (T ) that are the sets containing the edges forming the three portions of the boundary D , N and T . We define h K and h E to be the diameter of the element K and the length of the edge E respectively and we also define h max as the maximum of h K on the mesh T . Now we introduce the polynomial degrees for the approximation in our DG method. Hence, for each element K of the mesh T we associate a polynomial degree p K ≥ 1 and we introduce the degree vector p = { p K : K ∈ T } and we define p min as the minimum of p K on the mesh T . We assume that p is of bounded local variation on all meshes in the sense that for any pair of neighbouring elements K, K ∈ T , we have where ϱ ≥ 1 is a constant independent of the particular mesh. For any E ∈ E (T ) , we introduce the edge polynomial degree Hence, for a given partition T of and a degree vector p on T , we define the hp -version DG finite element space by where P p K (K) is the space of polynomials of degree at most p K and R is the set of rigid motions, see Remark 1 . The exclusion of rigid motions from the DG space is useful for the analysis. In practice this constraint is imposed as explained in Remark 1 .
Given an edge E ∈ E int (T ) shared by two elements K + and K − , we define denoting with + and − the values from the two elements. Moreover, we define n + K = (n + x , n + y ) the outward unit normal on the boundary ∂K + of an element K , then we define the jump [[ · ]] operator and the average { · } operator on vectors and tensors as: Thus, the DG approximation problem (1) reads as follows: where the bilinear form where · 0, K and · 0, E are respectively the L 2 -norm on an element K and on an edge E .

A priori convergence results
In this section, we prove the a priori convergence of the DG method. For this purpose we need to introduce an additional function space: Remark 2. From now on the notation a ࣠ b is used to denote a ≤ c b , where c is a constant that may depend on the coefficients in (1) , the value of γ and the constant in Theorem 2.1 . The constant c is always independent of the sizes of the elements and the orders of polynomials in the elements. Lemma 1. Let K ∈ T be a triangular element and u a function in [ H s ( K )] 2 , with s ≥ 1 . There exists a positive constant C π depending on s and on the shape regularity of the mesh but independent of u , p K and h K , and a polynomial π u ∈ [ P p K (K)] 2 , with where μ = min (p K + 1 , s ) and where E ⊂ ∂K.
This result for the scalar case is presented in Lemma A.7 in [14] . The definition of the projection operator is in [15] for the p -case and in [16] for the hp -case. The extension of the operator from the scalar case to the vector case is trivial and it consists in applying the operator to each component of the vector function u .
where h max and p min are defined in Section 4 and where μ = min (p min + 1 , s ) and p min ≥ 1 .
Proof. The proof is based on [14,Theorem 4.1] , but in order to apply it to linear elasticity, we extended it in several ways including to remove rigid motions. First we need to introduce a second DG norm: it is straightforward to see that ||| · ||| T ≤ · DG , T . Also Theorem 3.3 in [14] can be used to prove the continuity result in our case, i.e.
Moreover, Theorem 3.5 in [14] can be used to prove the coercivity result in our case, i.e.
v h Using the interpolation operator defined in Lemma 1 , we have: then by (17) and the orthogonality of the DG solution, Furthermore applying (16) we obtain The result comes directly from (18) noticing that The regularity assumption in Theorem 5.1 may be weakened using the same argument for the pure diffusion case as in [17] .

Reliability of the error estimator
This section contains one of the main results of this paper which is the proof of reliability for the proposed error estimator (19) . Theorem 6.1 guarantees that, up to a constant independent of the sizes of the elements or their order, the error estimator is an upper bound for the true error using the DG norm. This result is very useful in practice since the error estimator is computable also when the exact solution of a problem is not known, therefore it can be used to determine how accurate is the computed solution. Also, since the error estimator is computable, it can be used to drive adaptivity to improve the accuracy of the computed solution. In Section 8 a series of examples are presented to show the usage of the error estimator.
The error estimator presented in this paper is given by where the three terms under the sum are defined as where f h and g N , h are the L 2 projections of f and g N onto the finite element space and where g D , h and g T , h are approximated as in [18] of traces of functions in H 1 ( ) such that and such that g D , h , g T , h are on each edge E on the corresponding portions of the boundary the best approximations of g D , g T in the interpolation space [( L 2 ( E ), H 1 ( E )) 1/2 ] 2 ; we refer the reader to [19] and the references cited therein for details on the construction of the approximations and to [20] for the definition of the interpolation Sobolev space.
Due to the fact that data like f or the values of the boundary conditions may not be represented exactly in the finite element space, the analysis include an oscillation term. This term is high order compared to η err and for this reason is rarely computed in practice. Also the oscillation term is the sum of different terms: where Next, we introduce the reliability result, which is the main result of this section. In the repository [21] we included a document containing a longer and more detailed version of the proof that is too long to include here.
where C is a positive constant independent of the mesh nor the order of the elements used.
The rest of the section is devoted to the proof of Theorem 6.1 . The proof is based on [18] where a reliability proof for an error estimator for the Laplace equation is presented. However, to prove Theorem 6.1 for linear elasticity, the approach has been changed to deal with the fact that the model problem is a system of equations. Also, the presence of many different boundary conditions increases considerably the number of terms to consider in the proof.
In order to carry out the analysis we need to introduce an auxiliary continuous problem similar to (1) which is only used in the analysis and never in computations: Since linear elasticity is a linear problem and its solutions depended continuously on their data, we The introduction of this auxiliary problem is essential in the analysis to isolate the oscillation term. For DG methods, the construction of the upper bound for the DG norm using the error is done in two steps bounding separately the conforming and the non conforming part of the error. To this end, we define the space V c 2 which is a conforming version of the DG space. Then, we decompose the discontinuous Galerkin (13) . We also define V c p , 0 (T ) which is the subspace of V c p (T ) containing functions satisfying the following boundary conditions imposed strongly: u = 0 on D , σ(u ) · n = 0 on N , u · n = 0 on T and t (u ) · n = 0 on T . We also assume that there is an interpolation operator I hp : V p (T ) → V c p (T ) that satisfy the following inequalities: with E ⊂ ∂K . Examples of similar interpolation operators are the Scott-Zhang presented in [22] and used in [18] or the operator in [8] .
We can then split the solution as: with u c h ∈ V c p , 0 (T ) and u r h ∈ V ⊥ p (T ) , then using the triangle inequality and (26) , we obtain The first term on the rhs of (27) can be bounded using (22) and (13) noticing that u −˜ u is zero on the internal edges: Then, in order to obtain the upper bound for the conforming part of the error ||| ˜ = 0 in the interior of the mesh and so denoting we obtain from Theorem 2.1 : we have The term D ( ˜ u − u , v ) can be bounded using the Cauchy-Schwarz inequality and (28) by The next step is to bound the H 1 -seminorm of u r h using the jump over the faces. Similar results have been already used in other works, like in [18] . However, such results are for problems with only Dirichlet type boundary conditions, which means that they are not suitable in this context. We need to bound the H 1 -seminorm of u r h using only the jump over interior faces and faces on the portion of the boundary where only Dirichlet type of boundary conditions are imposed. Such result is proved for scalar problems in Theorem 2.1(ii) in [23] and for the case D ∪ T = ∅ , which is admissible for our model problem, we use Theorem 2.1(iii) from [23] . Applying the results from [23] to our problem, we have that there exists a projection operator π hp : where the faces on the boundary considered on the rhs are also the ones on T for just the normal component since for that component the boundary condition is of Dirichlet type.
Then, since any function v r h ∈ V ⊥ p (T ) can be seen as v r h = v h − π hp v h for some v h ∈ V p (T ) , we can apply (35) on the first term of (34) and have Also the other two terms can be bounded straightforwardly by K∈T η 2 J,K .
Then the term D (u r h , v ) can be bounded using the Cauchy-Schwarz inequality and Lemma 2 by To bound the remaining term D (u − u h , v ) we use the fact that problem (12) can be rewritten as Then, we have choosing v as in (31) and using (1) : The next lemma is fundamental to bound the conforming part of the error with the error estimator.

Lemma 3.
Considering u h the DG solution of problem (12) and for any continuous function v with v h := I hp v , we have:

Proof. Applying integration by parts to the first term in
The term K∈T ∂K σ(u h ) · n K · (v − v h ) ds can be further treated in the standard way for DG: where in the last step we used the fact that [[ v ]] = 0 in the interior of the mesh. Adding also the term K ( u h , v h ) we have: In view of this, the rhs in Eq. (38) can be split in four terms defined as: The rest of the proof consists in bounding each term separately. To bound T 1 , we use the Cauchy-Schwarz inequality, (23) and v h := I hp v : To bound T 2 , we use the Cauchy-Schwarz inequality, (25) , the shape regularity assumption on the mesh T and the facts that v = 0 on D and v · n = 0 on T : To bound T 3 in the interior of the mesh, we use the Cauchy-Schwarz inequality, (25) and the shape regularity assumption on the mesh T : In a similar way T 3 is bounded on the Neumann portion of boundary: On the Dirichlet portion of boundary, the term T 3 is null since v = 0 on D : On the traction portion of boundary we use v · n = 0 to bound T 3 : To bound T 4 in the interior of the mesh, we use (4) , the Cauchy-Schwarz inequality, the standard hp -version of the trace inequality [24] and the shape regularity assumption on the mesh T : where in the last step we used (24) : In a similar way T 4 is bounded on the Dirichlet portion of boundary: In a similar way T 4 is bounded on the traction portion of boundary: The statement of the lemma is a consequence of all the above bounds.
Finally the proof of Theorem 6.1 is achieved by constructing an upper bound of (27) using (28) , Lemma 2 and (39) .

Efficiency of the error estimator
In this section we prove the efficiency of the error estimator exploiting the properties of bubble functions as in [8,[25][26][27] . The efficiency result cannot be shown uniformly in the polynomial degree since inverse estimates optimal in the polynomial degree are not currently available.
The efficiency result states that it is possible to construct an upper bound for the error estimator using the error in the DG norm. This, together with the reliability result, establishes that the error in the DG norm and the error estimator are linearly proportional quantities up to the oscillation terms that are considered of higher order. In other words we have that the error estimator cannot be "too far away" from the real value of the error.
For an element K the bubble function ψ K is positive real valued function with compact support contained in K and bounded by 1 in the L ∞ norm. Similarly, for an edge E we define the bubble function ψ E as a positive real valued function with compact support contained in E , where E is the union of the elements with an intersection with E of dimension 1, and bounded by 1 in the L ∞ norm. In [27] the explicit definitions of ψ K and ψ E for the 2D case can be found: denoting by λ K , i with i = 1 , 2 , 3 the barycentric coordinate of the vertex i of K Denoting with i = 1 , 2 the vertices of E and assuming they are not hanging nodes and with λ K + ,i and λ K − ,i the corresponding barycentric coordinates for the elements K + and K − forming E , For edges touching hanging nodes, the construction using an auxiliary mesh presented in [28] can be used. The construction is very technical and for brevity it is not reported in here.

Lemma 4.
Bubble functions can be constructed such that the following results hold for any element K ∈ T , for any edge E ∈ E (T ) and where C is a positive constant independent of the mesh nor the order of the elements used.
Proof. Starting with η J , K , we have that [[ u ]] = 0 on all interior faces, then where ||| · ||| ω K is the DG norm computed on all the elements intersecting K and its faces. Similarly u = g D on D and u · n = g T · n on T , so we have Moving on to the term η R , K and assuming that it is non-zero, we define w := Then applying integration by parts and using the fact that w | ∂K = 0 we have Then using (4) , (41) and the fact that ψ K ࣠ 1, we obtain Finally dividing both sides by In case that η R,K = 0 , then any non-negative quantity can be used to bound it and in particular (45) In case that h E p E ( σ(u h ) · n − g N,h ) E, 0 = 0 , then any non-negative quantity can be used to bound it and in particular (48) holds. Finally to bound the term on T we define w := h E p E t (u h ) · n ψ E and we proceed in the same way as before using the fact that t (u ) · n = 0 and applying integration by parts. This leads to the bound Putting together the bounds for all terms, the proof of the theorem is concluded.

Numerical examples
In this section the a posteriori error estimator, (19) , for smooth and non-smooth problems will be shown to be efficient, reliable, and with an exponential performance in the error estimate value (19) and the error in the DG norm (13) .
The finite elements used in the simulations are arbitrary high order triangular elements as defined in [29] . The hpadaptive strategy used here was originally proposed in [30] and was shown to be proficient for finite elements in [31] . To uniformly refine in h set δ 2 = δ 1 = 0 and to uniformly refine in p set δ 2 = 1 and δ 1 = 0 . The strategy assumes that if η 2 K > δ 2 η 2 max , where η 2 max = max K∈T η 2 K , the real solution in element K is non-smooth and so h-refinement is necessary.
If δ 2 η 2 max ≥ η 2 K ≥ δ 1 η 2 max the real solution is assumed smooth however the polynomial order is too low to capture the real solution to a sufficient accuracy. Many adaptive strategies exist, a review is provided in [30] , this strategy was chosen as it was the simplest to facilitate and demonstrate the efficacy of the error estimator. To perform h adaptivity only, set δ 2 = δ 1 = 0 and to adaptively refine in p set δ 2 = 1 and δ 1 = 0. To adaptively refine in hp simply ensure δ 2 = δ 1 . All adaptive strategies are halted once the number of degrees of freedom (ndof) of the linear system exceeds 10 4 . Algorithm 1. hp -refinement strategy: for parameters δ 1 , δ 2 with 1 ≥ δ 2 ≥ δ 1 ≥ 0 .
2) Identify the set of elements to refine in p: 3) Increase p K by one for K ∈ T p _ re f .

4) Identify the set of elements to refine in
4) Identify any elements K ∈ T ∩ T , where T is the refined mesh, that will have more than one hanging node on a face and add to T hre f . 5) h refine all elements K ∈ T h _ re f to create the new mesh T . 7) Ensure for every pair of neighbours K, K ∈ T that (8) is satisfied, otherwise add K , where p K < p K , to the set T p _ re f . 9) Increase p K by one for K ∈ T p _ re f .

Smooth solution problem
We consider a small strain linear elastic problem on (x, y ) ∈ = (0 , 1) 2 , where x and y are in metres. The problem acts in plane strain with Young's modulus E Y = 5 2 Pa, and a Poisson's ratio ν = 1 4 . The right-hand side f of problem (1) is chosen such that the exact analytical solution is smooth over the entire domain, u = sin (2 π x ) sin (2 π y ) sin (2 π x ) sin (2 π y ) . The initial mesh is conforming and is constructed from 32 elements with p K = 3 for all K ∈ T . The right-hand side is applied to the problem as a body force and with g D = 0 on D where ∂ = D . The hp -adaptive strategy uses δ 2 = 0 . 7 and Fig. 1. Square domain -performance of the error estimator, ηerr (19) , and the error in the DG norm ||| u − u h ||| T against √ ndof for different adaptive strategies. δ 1 = 0 . 07 , from [31] . As noted in [31] , since the solution is regular over the entire domain, and therefore smooth, adaptive p -refinement would produce the greatest reduction in error per unit cost in degrees of freedom. However it is demonstrated here for linear elasticity that this hp -adaptive strategy is still capable of producing exponential convergence. Additionally for comparison, the h -adaptive strategy uses δ 2 = δ 1 = 0 . 07 . In Fig. 1 the error estimate value and the error in the DG norm are plotted against √ ndof for, hp -adaptive, h -adaptive, and uniform h -refinement on a linear log scale. Fig. 1 shows only the hp -adaptive strategy to achieve exponential convergence for the DG norm and error estimate value, this is demonstrated by the (roughly) straight lines. Additionally in Fig. 2 we plot ||| u −u h ||| ηerr , against refinement step for the hp -refinement strategy. The variation in is oscillatory within a small range, which supports the fact that η err is efficient and reliable for smooth problems.
The hp -strategy employed here will always perform some h -refinement, unless δ 2 = 1 . Other hp -adaptive methods can achieve exponential convergence though p -adaptivity only, these adapt by evaluating whether the solution is locally smooth on an element by examining the decay of the element's Legendre coefficients [32,33] . A thorough investigation is presented in [31] .

L-shaped non-smooth problem
The next problem considered is a linear elastic problem on (x, y ) ∈ = (−0 . 5 , 0 .   (19) , and the error in the DG norm, ||| u − u h ||| T , against 3 √ ndof for different adaptive strategies. is singular at (x, y ) = (0 , 0) , On the boundary of the problem ∂ = D we have g D = u on D . The initial mesh is conforming and is constructed from 6 elements with p K = 3 for all K ∈ T . For this problem the solution at the convex corner is non-smooth, the error estimator here is therefore likely to be higher here than in the remainder domain. As in [31] we set δ 2 = 0 . 7 to capture the regions elements where the solution is non-smooth and perform h -refinements. δ 1 = 0 . 07 is set to capture remain regions of the domain where the solution is smooth but not sufficiently refined. The h -adaptive strategy used δ 2 = δ 1 = 0 . 07 . On Fig. 3 3 , see [34] . For the singular problem the hp -adaptive strategy achieves exponential convergence of the error estimator and the error in the DG norm, this is demonstrated by the roughly straight line on the linear-log plot. Additionally, Fig. 4 shows the hp -strategy to refine in h around the singularity and p in regions where the solution is smooth, consistent with [31,35] . Last the oscillations in Fig. 5 , show the error estimator to be efficient and reliable for singular problems.

Crack in a plate problem
Last we consider a problem with a stronger singularity, a crack in a plate acting in plane strain with E Y = 5 2 Pa and and p r = 1 Pa is an applied pressure. This is a mixed mode crack problem with jumps in both components of u across the crack faces, represented by the line ([0, 0.5] × {0}) in metres [36] . The near crack tip displacement field was first derived by Irwin [36] . It is presented here in polar coordinates ( r , θ ), where the crack tip is the origin, and κ = (3 − 4 ν) , r a << 1 , with a = 0 . 5 m as length of the crack and K I and K II as stress intensity factors which are depending on the loading, geometry and boundary conditions of the problem. A displacement solution does not exist for the entire domain, however the stress singularity at the crack tip is clearly stronger than that found in Section 8.2 .
To model the problem, the initial mesh is conforming and is constructed from 6 elements with p K = 3 for all K ∈ T .
The hp -adaptive strategy used δ 2 = 0 . 7 and δ 1 = 0 . 2 , and the h -adaptive strategy used δ 2 = δ 1 = 0 . 2 . Although δ 2 = 0 . 07  produced exponential convergence with hp -adaptivity, δ 2 was raised to 0.2 as this produced faster convergence by preventing unnecessary p -refinement in areas of the domain where the solution could be relatively well represented by low order polynomial functions. In Fig. 6 the error estimator values are plotted against 3 √ ndof . The element and polynomial order distribution is presented in Fig. 7 . The highest h -refinement levels are generated at the singular stress field at the crack tip. hp -refinement also occurs in the top-left corner of the problem where the stress field is also non-smooth, this results in relatively high error estimate values however not as high as the crack tip. For this singular problem only the hp -adaptive refinement scheme produced exponential convergence, the h -adaptive and h -uniform schemes produced only polynomial convergence.

Conclusions
The paper has presented, for the first time, an hp a posteriori error estimator for the symmetric interior penalty discontinuous Galerkin finite element method for linear elastic analysis. The error estimator was incorporated into an hp -adaptive finite element solver and verified against smooth and non-smooth problems with closed-form analytical solutions as well as being demonstrated on a non-smooth problem with complex boundary conditions. The hp -adaptive finite element analyses achieve exponential rates of convergence and are contrasted against uniform and adaptive h refinement.
The paper has provided a complete framework for adaptivity in the symmetric interior penalty discontinuous Galerkin finite element method for linear elastic analysis. This will allow engineers and scientists to use the method to obtain highlyaccurate but efficient stress analysis results in areas where the displacement/stress solution is of paramount importance, fatigue analysis for example.