Comparative Analysis of C- and C-GFEM Applied to Two-dimensional Problems of Confined Plasticity

For many practical applications in engineering, a complex structure shows linear elastic behavior over almost all its extension, but exhibits confined plasticity contained in some small critical regions, e.g. stress concentrations in fillets and sharp internal corners. The behavior of C 0 and C k -GFEM is investigated in this class of problems. The first goal of this study is to verify the actual formulation of the C k -GFEM for two-dimensional elastoplasticity, as a modification of the C 0 -GFEM formulation. The C k -GFEM is based on a set of basis functions with C k continuity over the domain. The approximation functions are constructed from a C k continuous partition of unity, over which polynomial enrichment functions (or any special function) can be applied, in the same fashion as in the usual C 0 -GFEM. In this way, the finite element approximations show continuous responses for both displacements and stresses across inter-element interfaces. An investigation is performed to assess the behavior of higherregularity partitions of unity against conventional C 0 counterparts. The irreversible response and hardening effects of the material is represented by the rate independent J2 plasticity theory with linear isotropic hardening of material and von Mises yield criteria, being considered only monotonic loading and the kinematics of small displacements and small deformations. The focus herein is to enlighten any possible advantage of smoothness in the presence of plastification phenomena, seeking for improvements in capturing the evolution of the process zone.


Abstract
For many practical applications in engineering, a complex structure shows linear elastic behavior over almost all its extension, but exhibits confined plasticity contained in some small critical regions, e.g. stress concentrations in fillets and sharp internal corners. The behavior of C 0 -and C k -GFEM is investigated in this class of problems.
The first goal of this study is to verify the actual formulation of the C k -GFEM for two-dimensional elastoplasticity, as a modification of the C 0 -GFEM formulation. The C k -GFEM is based on a set of basis functions with C k continuity over the domain. The approximation functions are constructed from a C k continuous partition of unity, over which polynomial enrichment functions (or any special function) can be applied, in the same fashion as in the usual C 0 -GFEM. In this way, the finite element approximations show continuous responses for both displacements and stresses across inter-element interfaces. An investigation is performed to assess the behavior of higherregularity partitions of unity against conventional C 0 counterparts. The irreversible response and hardening effects of the material is represented by the rate independent J 2 plasticity theory with linear isotropic hardening of material and von Mises yield criteria, being considered only monotonic loading and the kinematics of small displacements and small deformations. The focus herein is to enlighten any possible advantage of smoothness in the presence of plastification phenomena, seeking for improvements in capturing the evolution of the process zone.
Comparative Analysis of C k -and C 0 -GFEM Applied to Two-dimensional Problems of Confined Plasticity

INTRODUCTION
Certain local characteristics of boundary value problems, such as high gradient, singularities and discontinuities, can be successfully modeled with the use of the generalized finite element method (GFEM), since it allows using a priori knowledge about the solution of a problem in the form of enrichment functions. Several studies have shown that GFEM has been used successfully in linear elastic fracture mechanics (Areias and Belytschko, 2005;Belytschko, 2001;Duarte et al., 2007;Laborde et al., 2005;Pereira et. al., 2010;Stazi et al., 2003;Tabarraei and Sukumar, 2008). However, a real structure is a very complex body with stress states whose values based on the linear elasticity may exceed the elastic limits. The works of Kim et al. (2009) and Kim et al. (2012) proposed a GFEM with global-local enrichment functions (GFEM gl ) to numerically generate proper enrichment functions for three-dimensional problems with confined plasticity. According to Duarte and Kim (2008), the GFEM gl can be potentially more efficient than the standard GFEM or FEM. Kumar et al.(2014a) proposed an adaptive coupled finite element (FE) and element free Galerkin (EFG) approach used to simulate the nonlinear behavior of materials from two-dimensional stable crack growth problems, where EFGM has been used only in a region near the crack, while FEM has been used away from the crack. In the work of Kumar et al.(2014b), XFEM has been extended to simulate stable crack growth by J-R criterion, using finite strain plasticity under plane stress condition. They concluded that XFEM is an effective tool for modeling stable crack growth in homogeneous and bi-materials as it does not require the conformal meshing.
The piecewise polynomial partition of unity (PoU) functions used in the conventional instance of the method may be not the most appropriate for some kinds of problems, as in case of singular enrichments (Torres et al., 2014). In this context, the C k -GFEM, which is quite similar to C 0 -GFEM, presents the high regularity of the approximation as an attractive feature. This smooth version of GFEM is a mesh-based discretization methodology which brings together higher regularity, flat-top property, globally defined coordinates, compact support and no restriction regarding the shape of elements and/or clouds.
The importance of C k -GFEM comes from the fact that it is more robust than other kinds of higher regularity mesh-based approximations. Some C 1 approximants applied to problems that require such regularity, as Hermitian functions commonly used in thin plate bending problems, for instance, are very restrictive and do not favor a systematic procedure for building such functions for general element shapes or to allow p-refinement. Some approaches that exploit some benefits of C k -GFEM approximations in plate bending problems were presented in Barcellos et al. (2009), Mendonça et al. (2011 and Mendonça et al. (2013).
In this context, the goal of this paper is to compare the C k -GFEM and C 0 -GFEM performances in modeling two-dimensional problems involving elastoplasticity, considering problems with stress concentration (e.g. L-shaped domain), i.e., situations where the plasticized zone is confined to one or a few regions of the body. These kinds of problems are important due the high gradient of deformation field which occurs along the boundary of the process zone, which is difficult to be represented with coarse meshes and conventional functions of FEM. Thus, the ability of the approximation functions with C k arbitrary inter-element continuity to build continuous stress fields can be useful.
The irreversible response and hardening effects of the material is represented by rate independent J 2 plasticity theory with linear isotropic hardening of material and von Mises yield criteria, being considered only monotonic loading and the kinematics of small displacements and small deformations.
The paper is structured as follows. Section 2 briefly describes the C 0 -GFEM, which is used as a starting point for the C k -GFEM formulation developed in section 3. Section 4 presents the key elements on the elastoplasticity model considered herein, section 5 describes the test problems and presents the numerical results. Section 6 develops several concluding remarks about the findings.

GENERALIZED FINITE ELEMENT METHOD (C 0 -GFEM)
The GFEM is a combination of the standard finite element method (FEM) with concepts and techniques typical of meshless methods. This method presents an aspect of nodal enrichment that, in many situations, reduces or avoids the need of localized mesh refinement, making it very attractive for several analyzes.
The GFEM relies on a mesh that is used to define a partition of unit (PoU) and a domain for the numerical integration over which the enrichment of the PoU functions is performed. The set of PoU functions is employed to ensure the inter-elementar continuity, providing conformity of approximations that are improved by nodal enrichment strategy.
The enrichment functions are linked to individual clouds (patches) of elements around an arbitrary node in order to improve the quality of approximation in that neighborhood. Thus, one has the possibility to enrich the approximation only in a region of the problem domain, due to the compact support of PoU, without mesh refinement (Duarte et al., 2000), (Barros et al., 2004). Moreover, the essential boundary conditions can be imposed exactly as in the standard FEM (Strouboulis et al., 2001).
To build the GFEM approximation functions it is considered, for example, a conventional mesh of finite elements defined by N nodes with Cartesian coordinates in the domain Ω. In plane problems, usually the elements are three node triangles or four nodes quadrilaterals (and analogously in three-dimensional meshes). If the enrichment is performed with relation to node , a generic cloud is defined as the union of all finite elements adjacent to this node. The set of the shape functions belonging to each element associated with the node , compose the function on the support of the cloud . The set of functions associated to all nodes on the domain form a partition of unity (PoU). The enrichment functions related to the node , are denoted by (with ) and represent a set of q + 1 linearly independent functions. Enrichment functions can be chosen as polynomial functions, harmonic, anisotropic or even functions that are part of the solution of the boundary value problem (Babuska et al., 2002;Stroubouliset al., 2000;Melenk and Babuska, 1996). The local approximation subspaces can be denoted as which may also be enriched according to an adaptive method. In this study, uniform polynomial enrichments are chosen such that , where stands for the space of polynomials of degree less than or equal to p. For example, for p = 2 (quadratic basis), the set of enrichment functions is where are nodal coordinates of an arbitrary node j and and are the cloud characteristic dimensions around the node in the directions x and y, respectively. Thus, the GFEM approximation functions associated to the node result from the extrinsic enrichment of PoU, i.e., multiplying the PoU function with support on the cloud by components of : This implies an increasing in the number of unknowns per node in relation with those associated with the PoU. The resulting approximation function contains features of both functions, that is, the compact support of PoU and the approximation feature of the enrichment functions . The structure of GFEM offers more freedom in the choice of approximation functions compared to the standard FEM. The extrinsic way of to incorporate the algebraic refinements enables to vary the enrichment along the domain without compromising the conformity, which is ensured by the compact support of the PoU.
The generalized global approximation for the displacement on , denoted as , can be written as a linear combination of approximation functions associated with each node. The component , for example, can be written as: ( 2) where is a vector containing the full set of approximation functions and is a vector containing the nodal parameters and associated, respectively, with the PoU functions and the enriched functions . The continuity of the function on the whole domain is guaranteed by the characteristics of the PoU used.

Model Problem for Two-dimensional Elasticity
Let a boundary value problem (BVP) defined in a linear elastic domain , where the strong form of equilibrium equations is given by (3) where is the vector containing the stress tensor components, is the vector of body es, and denote complementary parts of the boundary , where Dirichlet and Neumann conditions are defined, respectively; and are prescribed displacements and tractions, respectively, and is the unit outward normal to . The operator have their definitions adapted to the vector notation considered for (see Zienkiewicz and Taylor (2005) for example) and the superscript T indicates transpose. The Cauchy stress is related to the displacement by the linear elastic relation , where is the vector containing the strain tensor components end is the constitutive matrix of the material, which is positive definite.
The variational form of this problem can be stated as: Find such that: where and are Hilbert spaces of degree order 1 (standard Sobolev space of square integrable functions whose first derivatives are also square integrable in the Lebesgue sense) defined on the domain Ω. The variational operators are defined as: where is the vector of displacements, is the test function vector, is the vector containing the strain tensor components, and is the thickness of the elastic body (in the z reference direction) considered here as constant.
The Galerkin approximation in Eq.(4), similarly to FEM, results in: Find such that: , whit and are being subspaces of finite dimension, wherein the first space is generated by the approximation functions and the second one is of test functions; is composed according to Eq.(2) and the component , for example, of the is given by: where is a vector of nodal parameters similar to vector . The only difficulty in the GFEM implementation is that the enriched basis can be linearly dependent, so that the system of equations resulting from Eq. (7) is positive semi-definite. The linear dependence occurs when the PoU and enrichments are both polynomial functions. This problem can be avoided by careful choice of functions (Oden et al., 1998), constraints in PoU (Melenk and Babuska, 1996), or the system can be efficiently solved by the numerical strategy proposed in Duarte et al. (2000) and Strouboulis at al. (2001). k -and C 0 -GFEM Applied to Two-dimensional Problems of Confined Plasticity Latin American Journal of Solids and Structures 12 (2015) 861-882

ARBITRARY CONTINUOUS FUNCTIONS (C K -GFEM)
C k -GFEM is quite similar to C 0 -GFEM, except that it builds an arbitrarily smooth PoU over finite element meshes (Duarte et al., 2006). Although this is not a meshless method, it preserves several attractive features of meshless approximations, such as the high regularity of approximation and partition of unit property. The C k -GFEM is also important because it can be efficiently applied to higher-order problems, as the plate problems of Kirchhoff and Reddy, which require solutions of continuity at least C 1 . The PoU functions can also be built as C  functions in case of clouds with convex support and considering exponential edge functions (the original Edwards (1996) procedure). In clouds with non-convex support they can be up to k-times continuously differentiable at the concave nodes, with k arbitrarily large, and infinitely differentiable in the rest of the cloud.
The technique used for the construction of PoU functions in the problems under analysis (Section 5), is described below for clouds with convex support and exponential edge functions (Edwards method).
Let us consider a conventional triangular finite element mesh defined by N nodes with coordinates in the domain Ω, and a set of weighting functions , with j = 1, ..., N. Each one is associated to a cloud as its compact support. Using Shepard's formula (Shepard, 1968), each PoU function is obtained by: . (9) Herein, triangle elements are considered aiming to employ them as integration cells directly, noting that the procedure is quite general for different element shapes or cloud configurations.
It can be seen that the set is such that , and and any compact subset of Ω intercepts only a finite number of supports. Therefore, is a partition of unity and its regularity depends only on the regularity of the weighting functions constructed to ensure the required continuity.
The weighting functions with convex support can be built from the product of the cloud edge functions associated with the cloud and defined in terms of parametric coordinates normal to each edge n of de cloud (Barcellos et al., 2009 andDuarte et al., 2006) as follows: , where is the number of edge functions to the cloud . Thus, the Shepard equation is used to define a partition of unity from weighting functions continuously differentiable, using exponential edge functions. According to Edwards (1996), a cloud edge function is a strictly positive function inside the cloud which vanishes, together with all its derivatives, towards an edge, rendering a weighting function.
The edge function used in this work is exponential with unitary value at the cloud node and its decay rate is controlled by a parameter , where is the normal dis-tance from node j to the edge n. For these two conditions, was used the following function edge (Barcellos et al., 2009) where and is the distance between the position and the edge n, is the edge midpoint, and the normal vector to the edge directed into the cloud. and are positive arbitrary constants. However, numerical experiments have shown that the most appropriate values are and , as suggested by Duarte et al. (2006) and Torres (2012). Therefore, at cloud node the function has the value (15) which is the same for every edge n of the cloud and imposing the constraint , one gets . It should be remarked that different choices of PoU functions are possible (it depends on the choice of the edge functions (Barcellos et al., 2009)), leading to different kinds of approximation functions.
However, besides exponential functions, many functions meet the requirements of having derivatives null on the edge. For example, in Barcellos et al. (2009) polynomial edge functions were subjected to numerical investigations. In this case of polynomial cloud edge functions, differently, the continuity is limited independently of the cloud geometry.
Herein, exponential edge functions are used, following the procedure described in Barcellos et al. (2009) which guarantee C  () approximation for meshes with only convex clouds. On the contrary, in the case of non-convex clouds, it is necessary to replace the edge functions of a pair of non-convex edges with a single new edge function obtained through a boolean product (Rvachev and Sheiko, 1995) before performing the weighting functions computation, as proposed by Duarte et al. (2006). Therefore, the resulting PoU is at least k-times continuously differentiable, and the resulting approximation functions of the product of Shepard PoU with the enrichment functions will have the same continuity, since the enrichments are at least also C k .

TWO-DIMENSIONAL PLASTICITY
This section presents the equations that govern the classical plasticity in the two-dimensional context, which can be found in Chen (1988), Simo andHughes (1998) andSouza Neto et al. (2008).
The formulation adopted is based on the classical rate independent J 2 flow theory for small strain. Its main features are: von Mises yield criteria, linear isotropic hardening of material, hypothesis of associativity to the hardening law and normality rule for plastic flow. The Newton-Raphson was the iterative and incremental scheme adopted.
Initially the isotropic hardening will be used, however, kinematic hardening and cohesion models can easily be introduced to the algorithm. The plastic flow is regarded as an irreversible process and is characterized in terms of the history of the following variables: strain tensor , plastic strain tensor and isotropic hardening internal variable , related to the evolution of plastic deformation.
Considering the additive decomposition of the total strain tensor , the isotropic linear elastic constitutive tensor , the deviatoric stress tensor s and the yield stress , the equations that govern the model are: 1) elastic stress-strain relationship: ; 2) von Mises yield criterion: , where e is the modulus plastic of isotropic linear hardening; 3) flow rule: , where ; 4) hardening law: ; 5) accumulated plastic deformation rate: ; 6) Kuhn-Tucker complementary conditions: , , ; 7) consistency condition: (if , where is the yield function rate.
The parameter is a non-negative function, called consistency parameter, which represents the plastic flow rate satisfies the Kuhn-Tucker and consistency conditions.

NUMERICAL RESULTS
In this section, two numerical experiments are conducted. The goal is to compare the C k -GFEM and C 0 -GFEM performances for two-dimensional elastoplastic problems. The elastoplasticity algorithm was introduced in a code developed for linear elastic fracture mechanics (Torres, 2012).
The numerical experiments were conducted considering regular and uniform triangular meshes, with simplex elements, for which only convex clouds occur, in such a way that the functions are C  continuous for the smooth version of GFEM.
The algebraic refinement consists of uniformly applied polynomial enrichment functions of degrees p = 1 to 4. The numerical integration on the elements was performed by the Wandzura's symmetric quadrature in triangles (Wandzura and Xiao, 2003) with 175 points. The same quadrature was used for the C k -GFEM and C 0 -GFEM independently of the p-enrichment applied. It is worth noting that this amount of integration points is used here only for sake of evaluation, to reduce the errors due to integration from discretization errors to negligible levels. For integrations in the Neumann boundary was used the Gauss-Legendre quadrature with 25 points, for each element edge.

Internally Pressurized Cylinder
The first example is a classical problem of a long thick-walled cylinder subjected to a gradually increasing internal pressure. The dimensions of the problem, the material parameters and the finite element mesh adopted are shown in Figure 1. The cylinder is discretized by 64 axisymmet-ric elements. The maximum pressure, P, used to this problem was 0.18 GPa. Nine uniform pseudo time steps were used for incremental analysis. Figure 2 shows the radial displacement at outer surface of the cylinder versus applied pressure computed by the C 0 -GFEM and C k -GFEM for b = 1, 2, 3 and 4, that represents the polynomial degree of reproducibility of the proposed GFEM approximation subspace. Figure 3 and Table 1 show the norm of the relative error of the radial displacement at outer face versus degrees of freedom obtained via C 0 -GFEM and C k -GFEM for the maximum pressure level (P = 0.18 GPa). On notice that b = p +1 to C 0 PoU and b = p to C k PoU, where p is the degree of the polynomial enrichment, as identified in Mendonça et al. (2011). The relative error of the radial displacement is computed by , where and are the exact (Hill, 1950) and approximated displacements, respectively. It can be observed that the C 0 -GFEM and C k -GFEM results are close to the analytical value and exhibit similar errors for the same values of b. As expected, the worst values of the approximations of displacement occur for b = 1. According to Figure 2, the approximation of the displacement for b =1 worsens at pressures above the yielding pressure onset, whose analytical value is 0.103 GPa (Souza Neto et al., 2008). This suggests that for this problem (choice of mesh, PoU, enrichment function) both methods do not approach well the plastic solution when b = 1.
The circumferential stress obtained at integration points for P = 0.1 GPa (elastic solution) and P = 0.

L-shaped Domain
The L-shaped domain is a classical problem of the elasticity theory for which an analytic solution is known (Szabo and Babuska, 2011). An elastoplastic version of this problem is considered here aiming (a) to verify the implementations for a general two-dimensional loading condition, (b) to investigate the effects of smooth discretization in the approximate stress fields, and (c) to compare the evolution of a process zone at the reentrant vertex neighborhood for the two kinds of regularity considered.  The dimensions, material parameters, loading and mesh (Mesh 1) used in this example are shown in Figure 6. The analysis is carried out assuming plane strain conditions. The domain of the problem is discretized by 96 elements. Ten uniform pseudo time steps were used for incremental analysis. The numerical results are compared with a overkill solution from ANSYS ® (Reference 1) for which was used a very refined regular and uniform triangular (3 nodes) mesh with 60 000 elements and 60 799 degrees of freedom. It should be remarked that for the L-shaped domain problem, considering elastic material, an infinity stress value occurs at the reentrant corner for the minimum amount of applied external force. This means that in the physical problem, in case of elastoplastic material, the yield criteriais reached once any external force is applied. Of course, for discretized problems, this fact is generally not true.  The strain energy, defined as , is used as a global convergence measure to verify the numerical implementation. Another mesh with 384 triangular elements (Mesh 2, see Figure 7), was used in this verification, where the blue region is defined to allow the computation of a local deformation energy.  Table 2 shows the strain energy calculated throughout the domain and Table 3 shows the correspondent value for the blue region of the meshes shown in figures 6 and 7, which approximately contains the plastic zone. One can see that the results obtained by the C 0 -GFEM (Mesh 2) is greater than Reference 1 values, for some degrees b, which would indicate the reference solution is still inappropriate. For this reason the numerical results are also compared with the overkill solution from ANSYS ® which is similar to the first one but with triangles of 6 nodes (Reference 2), totaling 241 599 degrees of freedom. It is easy to see the present GFEM solutions converge to this new reference solution (Reference 2) as the degree b increases. Although the results for C 0 are closer to the reference value as the degree b increases, in both situations the results of C k for b = 1 are better, indicating some higher ability of the smooth functions of to capture the plastification for lower degree approximations. Latter, the evolution of the plastic zone will be investigated trying to confirm this hypothesis. Table 4 lists the number of iterations in the Newton-Raphson scheme required in each loading step obtained by C 0 -GFEM and C k -GFEM (Mesh 1), for b = 1, 2, 3 and 4. The table indicates that the total numbers of iterations required for both approaches are almost the same. The loading level for which the plastic flow threshold is reached has number of iterations different from 1. It is possible to note (see Table 4) that for b = 2 and 3 the plastic flow starts at the same loading level for both conventional C 0 -GFEM and smooth GFEM. However, for the lower approximation degree (b = 1) the C k -GFEM is capable of detecting hardening two load steps before the conventional counterpart, thanks to the derivatives of the C k -PoU, which are not constant inside the elements. Thus, it can be seen that in case of lower degree approximations, even for coarse meshes, the C k -GFEM captures the plastic flow in a lower loading level if compared to the C 0 -GFEM. In case of elastic materials, it is known the appropriate patterns of mesh refine-D G ments which should be applied close to these points. However, for more general kinds of geometric details, or in case of elastoplastic materials, the definition of an optimal mesh-refinement could be a difficult task. From these facts, it seems to be interesting to use higher regularity functions around reentrances, or other geometrical details. (/10 -4      The respective numbers of degrees of freedom used in C 0 -GFEM and C k -GFEM analyses for different values of b are shown in Table 5. Table 6 shows values of displacements in the y direction at points B and C for the Mesh 1 (see Figure 6) and the relative errors obtained by the C k -GFEM and C 0 -GFEM for different degrees b. It can be noted that the convergence of both forms of GFEM when the degree b increases is clear, and it is seen that only for b = 1 the values obtained by C k -GFEM exhibits smaller errors than those obtained by C 0 -GFEM. These results suggest that the approximation subspace generated by the C k -GFEM is "less flexible" than that for the C 0 -GFEM since, for the same number of degrees of freedom and higher degrees of b, the error obtained by C k -GFEM increases. This observation agrees with the understanding that C k  C 0 .  Reference 2 -0.5473 1.0473 Table 6: Values of displacements in the y direction at points B and C ( Figure 6) and relative error (in parentheses). Figure 8 and 9 shows pointwise values of equivalent stress over the line DG in the meshes 1 and 2, respectively (see Figure 6 and 7), obtained by the C k -GFEM and C 0 -GFEM, for b = 1, 2, 3 and 4. Whereas that in the C 0 -GFEM was performed an average of nodal values, clearly C k -GFEM and C 0 -GFEM are not confronted with the same rigor. In general, the results obtained by C k -GFEM are closest to the reference solution (Reference 2) when compared with the results obtained by C 0 -GFEM. For b = 3 and 4, the right side of the curves for C k -GFEM is already good and the greatest difference occurs at the extreme left for both meshes. k -and C 0 -GFEM Applied to Two-dimensional Problems of Confined Plasticity Latin American Journal of Solids and Structures 12 (2015) 861-882 The figure shows that better results could be achieved with a refined mesh. For example, the use of a graded mesh would capture singularity at a lower loading level, and therefore the plastic flow. Additionally, an algebraic refinement consisting of non-uniformly applied polynomial enrichment functions, such as considering higher b degree towards the reentrance, would also be used, or even another kind of enrichment function. However, the costs are reduced using a coarse mesh and smaller values of b, and in this case the C k -GFEM presents better results compared to C 0 -GFEM. It can be noted that for Mesh 2 and b = 4 the difference in equivalent stress at the extreme left point of the line DG is still very large, as the exact solution is singular, due to the linear hardening, and only polynomial basis functions were used, without any special refinement toward the reentrant corner. Then, for verification purposes a geometric mesh (see definition in Szabo and Babuska, 2011) with 72 triangular elements (Mesh 3, see Figure 10), was used in this verification. Figure 11 shows pointwise values of equivalent stress over the line DG in Mesh 3, obtained by the C k -GFEM and C 0 -GFEM, for b = 1, 2, 3 and 4. In general, the results obtained by C k -GFEM and C 0 -GFEM are similar, and for b = 3 and 4, the curves agree very well, even at the extreme left point. Note that although the exact solution is infinite at point D, the reference solution is finite at this point, since it is numerically calculated and uses a base that belongs to a space of regular functions. Figure 12 shows the relative error of equivalent stress at D point, with respect to the reference solution, versus degrees of freedom obtained via C 0 -GFEM and C k -GFEM. It is observed that the convergence in geometric mesh is much faster than in other ones, as expected. Figure 9: Values of equivalent stress over the line DG to Mesh 2 (see Figure 7). Some influence of smoothness can also be seen from the geometric aspect of the process zones for different degrees b at the end of loading. Figure 13 shows distributions of integration points where the yield condition was reached for the complete loading. This figure compares the C k -GFEM and C 0 -GFEM performances considering the high gradient of deformation field that occurs in the plastic zone. The results suggest the ability of the approximation functions with arbitrary inter-element continuity to represent the plastic zone more properly than the C 0 , mainly for lower degree approximations, as early mentioned (see reference solution in the Figure 14).

CONCLUSIONS
The aim of this study was to verify the GFEM implementation for two-dimensional elastoplasticity, and after that, compare through numerical experiments the C k -GFEM and C 0 -GFEM performances in problems with confined plasticity based on J 2 plasticity theory. For the two problems analyzed such comparison was performed using local and global convergence measures.
The results presented for the pressurized cylinder problem show that the quality of the stress is almost the same for both C 0 -GFEM and C k -GFEM, for the same degree of the approximation b. That is expected considering that this problem is essentially one-dimensional.
On the other hand, the results presented for the L-Shaped problem suggest that the C k -GFEM furnishes better approximations for the plastic region boundary and the equivalent stress values, mainly in the neighborhood of the reentrance, when compared with C 0 -GFEM. Furthermore, for b = 1 the C k -GFEM presents better results when compared to those obtained by C 0 -GFEM. This observation comes from the fact that the C k PoU functions have non-constant derivatives inside the elements, differently from C 0 tent functions, which favor capturing much localized features even using coarse meshes. Then, as can be seen in Table 4 (for b = 1), since the smooth PoU functions enable the detection of hardening for smaller load if compared to the C 0 continuous counterpart, it can be argued that the C k PoU functions around critical points are more appropriate for determination of a load level with cause plastic flow in the structure.
Finally, these are the first exploratory results towards a larger work which aims to use the C k -GFEM for composing local problem of the Global-Local GFEM (GFEM gl ). One of the next steps of this research is to keep the process of checking of the L-shaped problem, where investigations will be conducted to reduce costs by applying smooth PoU functions only in the region of the recess, using coarse mesh and small b. Figure 14: Contour of the plasticized region obtained from the reference solution for the complete loading.