On the sparsity of linear systems of equations for a new stress basis applied to three-dimensional Hybrid-Trefftz stress finite elements

Hybrid-Trefftz finite elements have been applied to the analysis of several types of structures successfully. It is based on two different sets of approximations applied simultaneously: stresses in the domain and displacements on its boundary. This method presents very large linear systems of equations to be solved. To overcome this issue, most authors have been careful in the choice of the approximation fields in order to have highly sparse linear systems. The natural choice for the stress basis has been linearly independent, hierarchical and orthogonal polynomials which typically result in more than 90% of sparsity in 3-D finite elements. Functions derived from associated Legendre and Chebyshev orthogonal polynomials have been used with success for this purpose. In this work the non-orthogonal polynomials available in the Pascal pyramid are proposed to derive a harmonic and complete set of polynomial basis as an alternative to the above-cited functions. Numerical tests show this basis produces accurate results. No significant differences were found when comparing the sparsity of the linear system of equations for both functions.


INTRODUCTION
The hybrid-Trefftz stress element formulation presents itself as an alternative for the dominant conforming singlefield based displacement element in computational analysis after the pioneering work of Pian (1964), Pian and Tong (1969) and de Veubeke (1980), which has been thoroughly compiled along with other non-conventional methods by Freitas et al. (1999).This formulation, considering the linear isotropic case for simplicity, consists on the independent approximation of the stress field in the domain of the element and the displacement field on its boundary.The Papkovitch-Neuber solution of Navier equation is used to derive the stress approximation fields to satisfy the Trefftz constraint, i.e., the displacement in the domain must satisfy locally all field equations.This technique imposes the use of harmonic potential functions for generating stress solutions.
The hybrid-Trefftz method, in particular when applied to 3-D problems, generates very large linear systems when high order polynomial approximations are used.For instance, a 50-element mesh with a hierarchical polynomial stress approximation of order 10 can create a matrix size larger than 20,000 × 20,000, making the sparsity feature paramount in the computation of the solution of the resulting linear system.The natural choice for the stress basis has been orthogonal polynomials to display high sparsity indices (Freitas, 1998), which typically results in high level of sparsity, say, more than 90% in three-dimensional finite elements.
In this work the non-orthogonal polynomials available in the Pascal pyramid are proposed as an alternative to orthogonal polynomial bases to be used as stress approximation functions in the context of 3-D hybrid-Trefftz elements.Homogeneous Harmonic Polynomial (HHP) functions derived from the Pascal's pyramid of polynomials proposed by Wang (2002) are applied to the Papkovitch-Neuber solution of the Navier equation to derive a complete set of 3-D stress and displacement bases.This procedure was applied with success to the analysis of plates and shells with hybrid-Trefftz elements by Martins et al. (2018).
The sparsity levels of the finite element matrices produced by this approximation functions are compared to those produced by one set of orthogonal and harmonic functions.Without loss of generality, in this work the basis derived from associated Legendre and Chebyshev (LC) polynomials is the chosen one.Both Legendre and Chebyshev stress approximation bases are orthogonal in [-1, 1] domain.It has been used with success by Freitas and Bussamra (2000), whereas the drawback of this choice is that in 3-D formulation the completeness of this stress basis is limited to the sixth degree.In addition, the accuracy of the stress predictions of the proposed element is analyzed through numerical tests.
Hybrid-Trefftz elements with both nodal and generalized variables framework have shown good performance in linear elastic (Freitas and Bussamra, 2000) and elastoplastic (Bussamra et al., 2001) analysis of solids with LC functions.In crack analysis, singular stress fields and stress concentration problems were analyzed using LC and Airy functions (Bussamra et al., 2014) and Kaczmarczyk and Pearce (2009), respectively.In multisite cracked solids, Chebyshev functions were applied in a nodal framework (Argôlo and Proença, 2016).Some authors applied the hybrid-Trefftz formulation to problems other than linear elastic mechanics.Fu et al. (2011) analyzed heat conduction in functionally graded nonlinear anisotropic materials using a nodal hybrid-Trefftz element.Cao et al. (2013), Lee et al. (2010) and Souza and Proença (2009) used complex variables derived from works from Muskhelishvili (1953) and Qin and Wang (2008) to approximate the domain fields when analyzing micromechanics of heterogeneous composites, crack singularities and the effect of selective enriching approximation functions, respectively.Wang et al. (2014) applied the dual reciprocity method to orthotropic potentials modeled with hybrid-Trefftz finite elements, dividing the solution into homogeneous and particular parts.Petrolito (2004) used bi-harmonic polynomials as approximation functions implemented in triangular elements with hybrid-Trefftz formulation to analyze stability and buckling of thick and thin plates in a 2-D approach and later studied vibration and stability on thick orthotropic plates with complex conjugate harmonic polynomials (Petrolito, 2014).Karkon and Rezaiee-Pajand (2016) studied thick orthotropic plates, with the difference of using orthotropic Timoshenko beam interpolation functions for approximation at the boundary fields and using both triangular and rectangular hybrid-Trefftz elements in various benchmark tests from the literature.Karkon (2015) also proposed triangular and rectangular hybrid-Trefftz elements to analyze anisotropic laminated plates.
No study was found, to the best of the authors knowledge, that compares not only numerical results but also the proposed function's numerical applicability measured in terms of the sparsity of their linear system of equations.Therefore, this work proposes to compare the aforementioned harmonic function sets in terms of sparsity and completeness.and Structures, 2020, 17(7), e307 3/17

FORMULATION
The hybrid-Trefftz stress element formulation derived in this work is based on the linear elastic fundamental governing equations, applied to a system with domain V and enclosed by a boundary Ć, referred to a Cartesian coordinate system: where vector ó and å gather the independent components of the stress and strain tensors in the equilibrium and compatibility equations, Eqs. ( 1) and ( 2) respectively; b � represents the prescribed body forces vector and u the displacements vector.The constitutive equation Eq. ( 3) is represented in the flexibility format with f being the flexibility matrix, symmetric and with constant entries when a linear, reciprocal elastic law is assumed.ó � č and å � č represent the residual stress and strain vectors, respectively.For simplicity, ó  , å  and b � are set to zero.Equation (4) stands for the Neumann boundary condition, applied in the static section of the boundary (Ć ó ), where t ̅ Ć are the prescribed tractions.Equation ( 5) stands for the Dirichlet boundary condition, applied to the kinematic portion of the boundary (Ć u ), where the displacements (u � Ć ) are prescribed.D is the differential equilibrium operator and D * is its Hermitian transpose.Both are linear and adjoint in the context of geometrically linear models.Matrix N contains the unit outward normal vector associated with the operator D.

Approximation fields
The element formulation used in this work is the hybrid-Trefftz.The term hybrid means two independent field approximations are made.One field is approximated in the domain, and the other on its boundary.Since the element formulated is of the stress model, the generalized stresses are directly approximated in the domain and the displacements in the boundary.The stress approximation is while the boundary displacement approximation is where S and Z contain the approximation stresses and displacements and X and q their unknown weights, respectively.

Trefftz constraint
The Trefftz constraint is enforced in the domain equilibrium stress approximation Eq. ( 6) by requiring it to satisfy locally the system of differential equations.This requirement renders the following condition: which means that S must represent a self-equilibrated stress field.Equation (1) can also be written in terms of the domain displacement, generating the well known Navier's equation of compatibility.It is obtained by substituting the compatibility equation Eq. ( 2) and the constitutive relation Eq. (3) written in terms of rigidity into the equilibrium equation Eq. ( 1), following: Latin American Journal of Solids and Structures, 2020, 17(7), e307 4/17  = kD T u in V, and (9) The Trefftz constraint is based on a self-equilibrated approximation field S directly associated to the domain displacements u.The displacements u in the domain is approximated by where U collects the functions associated with the displacement basis, X is the displacement vector and u Ć collects the rigid-body motion.Substituting Eq. ( 11) into Eq.( 9) results in the stress basis S, based in the domain displacement field

ON THE CHOICE OF U
According to Fu et al. (2012) some 3-D isotropic elasticity fundamental analytical solutions are available in the literature, as for instance the Boussinesq-Galerkin (Wang, 2002), Papkovitch-Neuber and quasi Hu Hai-Chang (Hu, 2008).However, Fu et al. (2012) reached the conclusion that from these presented solutions only a modified version of Papkovitch-Neuber's proposed by Wang et al. (2012) is able to directly formulate a linear independent and complete set of displacement approximation functions.Papkovitch (1932) and Neuber (1934) independently proposed a threedimensional solution to the Navier equation Eq. ( 10) for isotropic materials, which has the form where ϕ and Ψ are a scalar and vector harmonic displacement potentials, respectively, r is the position vector, ∇ is the gradient operator, ķ is the Poisson ratio and G is the shear modulus.Naghdi and Hsu (1961) and Mindlin (1936) shown that the Papkovitch-Neuber solution can provide complete solution of the Navier equation.This means the solutions are capable of representing every elastic displacement field possible in a three-dimensional problem.However, Eq. ( 13) can provide redundant solutions, and therefore not unique, so it is necessary to verify for linear dependencies.The choice resides in which harmonic potential is used to generate the displacement field to be substituted into Eq.( 12).Four independent set of functions exist in the proposed solution (Ψ 1 , Ψ 2 , Ψ 3 and ϕ), each of them to be substituted for the desired harmonic function set.
Moreover, Papkovitch (1932) also claimed the displacement potential ϕ could be set to zero without compromising the generality of the solution.Neuber (1934) claimed that any of the four harmonics could be set to zero with the same effect described by Papkovitch, but both statements were showed unsupported and inconclusive (Eubanks and Sternberg, 1956;Sokolnikoff, 1956), generating great discussion over the exact conditions that the four harmonic potentials could be turned into three.According to the investigation performed by Eubanks and Sternberg (1956) and followed by Naghdi and Hsu (1961) and Cong and Steven (1979) over the generality of the Papkovitch-Neuber potential, some of the conclusions found were: • if the domain is convex, any of the harmonic functions in Ų could be set equal to zero (regardless of the value of ķ) without loss of completeness; • the scalar function ϕ can be dropped if: o the domain is finite and star-shaped with respect to the origin; Eubanks and Sternberg ( 1956) demonstrated through a counter example that if 4ķ is an integer, ϕ cannot be dropped; • if 4ķ is non-integer and ϕ is a harmonic polynomial in x, y and z then ϕ can be dropped regardless of the form of the domain.Solids and Structures, 2020, 17(7), e307  5/17 The finite element geometry proposed to be used in is this work is a hexahedron, a convex three-dimensional domain.According to the conclusions presented above the potential ϕ can be dropped without any loss.Therefore, if potential ϕ is dropped the Papkovitch-Neuber solution Eq. ( 13) applied to Eq. ( 9) can be expressed by

R E T R
with: x 3 ) the point coordinates and ∂ 1 , ∂ 2 , ∂ 3 the partial differential operators in the three cartesian directions.

Legendre and Chebyshev harmonic potentials
Freitas and Bussamra (2000) assigned a polynomial stress basis derived from Legendre and Chebyshev hierarchical and orthogonal polynomials.The sparsity levels they are able to generate is paramount to the hybrid-Trefftz element due to the high number of degrees of freedom each element has, which leads to large linear system of equations.Their affinity with p-refinement is also an important and desired feature in hybrid formulations as it can exploit hierarchical function sets.Legendre and Chebyshev potentials generate complete sets of approximation functions, but only up to the sixth degree, as showed by Freitas and Bussamra (2000).The Legendre polynomials are given by The potentials generated from these functions are: It can be proven that this potential is harmonic if P n is a Legendre polynomial of degree n by performing Latin American Journal of Solids and Structures, 2020, 17(7), e307 6/17 The Chebyshev harmonic potentials were proposed by Freitas and Bussamra (2000), and can be separated into two subdivisions: Chebyshev ö and Chebyshev ϕ potentials, as follows: φ s y = r n sen(nθ) φ s x = r n sen(nθ) ϕ s z = zr n sen(nθ) ϕ s y = yr n sen(n) Legendre polynomials generate 3 different harmonic sets of functions producing a maximum of 9 independent fields.Chebyshev's generates 12 harmonic sets which produce a maximum of 36 independent fields.Together, they produce a maximum of 45 possible fields for each degree of approximation but there may exist linear dependent modes and these must be eliminated.This information is shown later in Tables 1 and 2.

Homogeneous Harmonic Polynomials
Aiming to build a harmonic polynomial set of independent functions derived from Pascal's pyramid trinomial distribution, Wang et al. (2012) applied the Laplace operator (the condition for a function to be harmonic) to the following polynomial, for a given degree of approximation n: As an example, consider n = 2, where a i , i = (1, …, 6) are constant coefficients, then f = a 1 x 2 + a 2 xy + a 3 xz + a 4 y 2 + a 5 yz + a 6 z 2 . (23) Applying the 3D Laplace operator: The result in Eq. ( 24) implies there is only one restriction, or dependency, among coefficients a i .Substituting Eq. (24) into Eq.( 23) it is possible to notice five independent terms arise.Factoring them according to the coefficients   , i = (1, …, 5) results in the five independent terms of this polynomial set, namely The main advantage of using this harmonic set of polynomials as approximation functions lies in its full completeness for every desired degree of approximation n, as shown by Wang et al. (2012) and Martins et al. (2018).

FINITE ELEMENT MATRICES
As stated by Freitas (1998) and Bussamra et al. (2001) there are different approaches to establish the finite element equations from the fundamental relations Eqs.(1-5) and the basic field approximations Eqs.(6-7), namely the duality, principle of virtual work and well-established variational statements.Here the virtual work approach is followed.
The element is based in the virtual work equation: This requires the stress field approximation ó to be in point-wise equilibrium within the element, and the boundary displacement field approximation u � to be the same along adjacent elements' common boundaries.As mentioned before, the generalized body forces are considered absent in this work.Substituting Neumann Eq. ( 4) and Dirichlet Eq. ( 5) conditions into Eq.( 24) follows: Taking the first variation in terms of the generalized stresses of Eq. ( 26) leads to Substituting approximations Eqs. ( 6) and ( 7) into Eq.( 27) follows As the solution is not trivial, δX ≠ 0. Substituting the constitutive relation Eq.
(3) into Eq.( 28) the following system of equations rises: Alternatively, the element equation Eq. ( 29) can be represented in the following form: where, Analyzing the equilibrium in the static boundary, given by the Neumann condition Eq. ( 4), and considering that the virtual work of the internal forces must be equal to the virtual work of the external forces   =   the following equation is obtained: Applying the domain stress approximation Eq. ( 6) and the boundary displacement approximation Eq. ( 7) in the following form Latin American Journal of Solids and Structures, 2020, 17(7), e307 8/17 into Eq.( 34) results in with Equations ( 30) and (35) render the system of equations in matrix form shown below: Matrix N is the normal operator related to the differential operator D. Vectors v and Q depend on the geometry of the element due to the prescribed displacements and prescribed tractions on its surface, respectively.Matrix A, vectors v and Q must be calculated for each unconstrained face.Matrix F of a given element is a square matrix of size equal to the number of accumulated stress fields generated by Papkovitch-Neuber solution.For each degree n Eq. 14 generates an approximation field that composes the stress approximation matrix S as shown in Eq. 38.The only exception is S 0 , which is a 6x6 identity matrix.

S = [S 𝟎𝟎 , S 1 , S 2 , …, S n ]
(38) Matrix Z contains the boundary displacement approximation basis defined in Eq. ( 7).Its functions are built from simple binomial distribution, defined in each face's local coordinates (ī 1 , ī 2 ).Considering n as the boundary displacement approximation degree, the number of fields obtained is given by Eq. ( 39), where each row represents the displacement approximation in one specific direction.
For n = 2: Equation ( 40) represents an unconstrained face.As an example, to simulate a simply supported case in a given face one of the rows representing the desired constraint direction should be removed, maintaining the remaining fields.In case of a clamped face, where movement in the three directions are restricted, the whole face is left out of the calculations and it is not accounted for in Eq. ( 31) nor in Eq. ( 36).This approximation does not satisfy face-to-face and edge continuity between elements.This effect is lessened as the exact solution is approached.An advantage of having this non-conformity is the higher continuity obtained in the stresses, which is desired in a stress element Freitas (1998).
In general, matrices F, A and vectors v, Q are defined for each element, and all of them depend either on the stress approximation S or displacement approximation Z.

COMPLETENESS OF THE STRESS BASIS
Observing Eq. (37), F is the only part of the linear system strictly dependent of the stress approximation S. Since the face integral A has influence from the boundary displacement Z, which is commonly constructed with simple independent polynomials, the sparsity analysis of integral A can be inconclusive in terms of the effects each proposed stress approximation functions have.
As h-and p-refinements are applied to increasingly complex problems, the linear system given by Eq. ( 37) can become very large and cumbersome to solve.Therefore, a common base for comparison must be established in order to Solids and Structures, 2020, 17(7), e307 9/17 obtain meaningful results regarding the sparsity levels of the LC and HHP functions.In order to do that, each of this function sets are evaluated in its completeness and in the independency of the generated approximation fields.

Completeness and independency of the stress fields
Considering that a three-dimensional complete field approximation of degree n, which has been defined by trinomial distribution in Eq. ( 22), has its size given by equation it is necessary to subtract from this group of possible solutions the ones that do not fulfill the imposed restrictions in Navier's equation The resulting number of accumulated independent fields after eliminating six rigid-body terms is given by Since a differentiation is applied at the Papkovitch-Neuber solution, the degree of S relates to the degree of n in the following manner Table 1 shows, for each degree n, the number of independent approximation fields that are obtainable from Eq. 14 and the accumulated amount of approximation fields carried over from the previous degrees, which is given by Eq. 41.The Papkovitch-Neuber solution provides 45 approximation fields for each degree of S. Out of the 45 obtained fields there exists linear dependencies among them which once eliminated will form the subset of linear independent fields that are used in Eq. ( 38).
From Tables 2 and 3, it is possible to note that the approximation basis formed by: • Legendre + Chebyshev ϕ is complete up to 1 st degree; • Legendre + Chebyshev φ is complete up to 3 rd degree; • Chebyshev φ + Chebyshev ϕ is complete up to 4 th degree; • Legendre + Chebyshev φ and ϕ is complete up to 6 th degree.
Table 2 shows how each orthogonal harmonic function set behaves when observed separately, in terms of independent approximation fields.In Table 3 these function sets are combined and have its completeness analyzed.It is possible to notice that the LC associated functions is able to generate complete stress approximation fields, according to Eq. ( 41), only up to the sixth degree.In the other hand, the HHP functions are able to generate complete sets of independent fields for any desired degree (Wang et al., 2012).

NUMERICAL IMPLEMENTATIONS
Each element's local coordinates are mapped on a master element, an 8-node hexahedron element, through a set of trilinear isoparametric functions.Its natural coordinates axes defined as [q, r, t], ranging from -1 to 1 from face to face and with its origin situated in the center of the element.The transformation functions are defined as where   ,   and   are the i th node coordinates of the master element.These functions are applied as shown below, where x local are the element's coordinates in each element's local system, and x master are the master's hexahedron element coordinates All surface and volume integrals are exactly calculated in the domain [1, -1] by using the Gauss-Legendre quadrature, as suggested by Zienkiewicz et al (2013).

SPARSITY RESULTS
Both LC and HPP potentials are compared in terms of the resulting linear systems of equations' sparsity.As previously mentioned, hybrid-Trefftz finite element analysis may generate very large linear systems, where the sparsity feature is more than desired.In this work sparsity is defined as the number of null terms in relation to the total number of terms present in the linear system.A 99% sparsity level implies that only 1% of the linear system terms are non-zeros.The following numerical examples were programmed using the software MATLAB© 2019b.Its finite element matrices are transformed to its sparse form through the sparse function, and the mldivide function is used as solver for the linear system of equations.

Bi-clamped beam under distributed bending load
A bi-clamped beam under a distributed bending load q = 1 applied on the top face is used to verify the accuracy of the results obtained by both proposed approximation basis, and to observe how the sparsity levels change when more Solids and Structures, 2020, 17(7), e307 11/17 elements are added.The beam has length of L and a square section of 0.2L, as shown in Fig. 1.The material is considered linear and isotropic with Young's Modulus E = 1 and Poisson's ratio ķ = 0.2.Three sets of degrees of approximations are used: [5,2], [7,3] and [9,4], where the first number is the domain stress approximation S degree and the second is the boundary displacement approximation Z degree.The displacement on the middle-bottom of the beam is evaluated and compared to the result obtained by a commercial finite element software analysis with 22,500 quadratic hexahedron displacement elements, v = 5.509qL/E.Table 4 displays the size and sparsity of matrix F (Eq. 31) of each element.This matrix is calculated only through the stress approximation function S, which results in a direct assessment of the sparsity level these functions can generate.In addition, the size of F is an important information since expressive computing time can be saved in the assembly of matrices F exploiting the fact that elements with the same material and size have the same F.In addition, the assembly the linear solving system (Eq.37) is well suited to parallel processing.These two proprieties were not exploited in this work.The size and sparsity of the resulting linear system of equations (Eq.37) is shown in Table 5.It is possible to verify that LC stress functions generated a slightly higher level of sparsity when compared to the HPP (the greatest difference, 1,42%, is found in the 1-element mesh of [9,4]

degrees).
Table 5 also presents the computation time spent to solve the linear system of equations.The increasing of the approximation degree (p-refinement) has more impact than h-refinement.It is also possible to note that HHP has, in most cases, a higher processing time.The two main reasons for that is the higher number of non-zero terms (lower Latin American Journal of Solids and Structures, 2020, 17(7), e307 13/17 sparsity) and since HHP is a complete basis for every approximation degree there are more equations in the resulting linear systems.This analysis ran with an Intel® Core™ i5-8400 CPU @ 2.80GHz processor.
By analyzing figure 2, it is possible to observe that p-refinement had a higher impact when less elements were present.When h-refinement was performed, very good results were obtained with few elements for the approximation sets used.Another remark is that both potentials reached very similar results.Taking LC approximation results as a benchmark, Table 6 shows the relative difference compared with HHP functions.

Cracked flat plate under traction load
In the subject of crack analysis, stress concentration and singular fields, the hybrid-Trefftz elements were applied by Bussamra et al. (2014Bussamra et al. ( , 2016) ) in a generalized framework with associated Legendre and Chebyshev polynomials, and by Kaczmarczyk and Pearce (2009) in a nodal framework with Airy functions.To verify the accuracy of the stress predictions of the proposed finite element, a structure with high level of stress gradient is analyzed.In this Section, a p-and h-refinement analyses of the results of stress intensity factor (K) of a cracked flat plate are shown.The plate is under a uniform far field tension ó = 1, with a crack oriented 90° from the stress application direction.The material is considered homogeneous and isotropic, with E = 1, ķ = 0.3, with dimensions according to Fig. 3.
A cracked plate behavior can be described by the stress intensity factor (K), which defines the crack tip stress state with help of the energy release rate (ÄG) that is related with the variation of deformation energy (ÄU) in the process of crack growth (Äa), say (Tada et al., 1973): Therefore, the stress intensity factor can be calculated through Eqs. ( 44) and ( 45) by increasing the crack length and calculating the variation of the energy release rate.The result can be compared with the one obtained by Tada et al. (1973), given by where F(a/H) for the given geometry is given by The analysis of the proposed beam was made using only the HHP functions.Approximation degrees of [5, 2], [7, 3], [9,4], [10, 3] and [10, 4] were used, where the first term is the domain stress approximation degree and the latter is the boundary displacement approximation degree.Five coarse meshes were analyzed, with 4, 12, 24, 48 and 72 elements (Fig. 3).For the given geometry and crack length the solution provided by Tada et al. (1973) is K = 5.0812.Results are displayed in Figure 4 and Table 7.

CONCLUSIONS
Orthogonal set of polynomials have been used to derive stress bases in hybrid-Trefftz finite elements by many authors.Legendre and Chebyshev (LC) potentials generate complete sets of orthogonal approximation functions, but only up to the sixth degree.In this work, a new stress basis (HHP) derived from the non-orthogonal polynomials from Pascal's pyramid is proposed for generating a 3-D hybrid-Trefftz finite elements.The numerical results obtained are compared to the results from Legendre and Chebyshev orthogonal and hierarchical stress approximation basis in terms of accuracy and linear system of equations sparsity.
The results showed that the HHP stress basis produces accurate stress and displacement approximations.When compared to the results produced by LC functions in the first example, the vertical displacement in the center of the bottom face of the beam varied less than 1,5% in maximum when compared with a benchmark solution obtained from displacement finite elements found in commercial softwares, but most results were around 0% when considering 4 decimal places.In terms of sparsity, HHP produced a very sparse linear system.It stayed slightly below the levels of the LC functions but the sparsity levels generated are close enough to justify its use, along with very good accuracy of the results.
The second numerical test showed that coarse meshes can produce very good stress intensity predictions.Approximation sets [5, 2], [7, 3] and [9, 4] produced good results, with less than 3% of error relative to the analytical solution and approximations [10, 3] and [10, 4] displayed the best results, with less than 1% of error when 24 elements or more are used.
In conclusion, HHP functions are a valid option to derive the stress approximation basis for the hybrid-Trefftz formulation, as it produces very sparse linear systems, showed good results, is complete for every approximation degree and its polynomial terms are hierarchical and easily generated.
On the sparsity of linear systems of equations for a new stress basis applied to three-dimensional Hybrid-Trefftz stress finite elements Felipe Alvarez Businaro et al.Latin American Journal of Solids

Figure 4 :
Figure 4: Convergence analysis of cracked plate with HHP functions.

Table 1 .
Expected number of independent fields of functions under the Trefftz constraint.Gray areas represent incomplete degrees.On the sparsity of linear systems of equations for a new stress basis applied to three-dimensional Hybrid-Trefftz stress finite elements Felipe Alvarez Businaro et al.

Table 2 .
Number of independent fields out of a total of 45.Each set of functions evaluated separately.Gray areas represent incomplete degrees.

Table 3 .
Number of independent fields out of a total of 45.Functions evaluated in pairs and in trio.Gray areas represent incomplete degree.
*Incomplete stress basis
*Incomplete stress basis

Table 6 .
Relative difference between LC and HHP approximation functions results.

Table 7 .
Stress intensity factor K for cracked plate using HHP stress approximation functions.