Thomas Elguedj

Yuri Bazilevs

Victor M Calo

Thomas J R Hughes

B and F Projection Methods for Nearly Incompressible Linear and Nonlinear Elasticity and Plasticity using Higher-order NURBS Elements

Keywords: Isogeometric Analysis, NURBS, Volumetric locking, Nonlinear Elasticity, Plasticity, Incompressibility, B-bar method, F-bar method

published or not. The documents may come

B-bar and F-bar projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements

1 Introduction

Preamble

Nonlinear finite element structural analysis is dominated by the use of loworder "displacement" elements that are specially designed to avoid volumetric or incompressible locking. Deformations that involve very small volume changes occur for rubber-like materials, and the elastic-plastic response of metals and undrained soils. Elastic-plastic analysis is ubiquitous in many engineering disciplines and, thus, it is essential that finite elements are able to accurately represent nearly incompressible deformations. Standard elements have difficulty or fail entirely in this situation. A full discussion of the early experiences and development of effective methodologies is presented in Hughes [START_REF] Hughes | Generalization of selective integration procedure to anisotropic and nonlinear media[END_REF], Chapter 4. It would seem at this stage of the development of finite element technology that higher-order approaches would play an important role in nonlinear structural mechanics, but this is not the case. The only approach that claims any success is the p-method, in which the polynomial order within elements is increased on a fixed mesh (see [START_REF] Szabó | Finite Element Analysis[END_REF][START_REF] Szabó | The p-version of the finite element method[END_REF]). The claim for standard higher-order C 0 -continuous finite elements is that volumetric locking is alleviated as the element polynomial order is increased; see, for example, Düster et al. [START_REF] Düster | p-fem applied to finite isotropic hyperelastic bodies[END_REF] and Heisserer et al. [START_REF] Heisserer | On volumetric locking-free behaviour of p-version finite elements under finite deformations[END_REF]. This seems to be the case, but there is evidence that the accuracy of the solution at any fixed polynomial order is far from optimal (some confirmation of this observation will be presented herein). In addition, numerical experience indicates standard higher-order elements are much more "fragile" than low-order elements. This lack of robustness is particularly apparent in nonlinear dynamic analysis of structures involving contact and impact, and subject to high wave number inputs, such as blast waves. The reasons for this have never been fully explained. However, recent vibration and wave propagation studies (see Cottrell et al. [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] and Hughes et al. [START_REF] Hughes | Discrete approximations in structural dynamics and wave propagation: p-fem vs. k -nurbs[END_REF]) have revealed that the higher modes produced by the p-method diverge with p. That is, whereas formal accuracy is increased, the improvement is confined to lower modes, while at the same time the higher modes get worse as p increases. This may explain why robustness decreases with p.

An alternative higher-order approach has recently emerged from Isogeometric Analysis (see Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]), namely k -refinement, in which discretizations of order p achieve C p-1 continuity on "patches" (roughly speaking, subdomains). It has been shown that k -refined meshes behave entirely differently than standard finite element methods with respect to higher modal components. In fact, in some cases all discrete modes converge to exact ones and nearly spectral accuracy is achieved. To us, this suggests that robust and higher-order accurate finite element methods applicable to nonlinear structural analysis may be a possibility, and this study represents a first-step in this direction.

At the same time, we need to deal with the locking problem in both small-and finite-deformation regimes. We feel the most widely applicable and practically useful approach is by way of a pure displacement formulation (i.e., no pressure degrees of freedom) as is generally employed in large-scale structural analysis programs (see for example [START_REF]LS-DYNA Keyword User's Manual[END_REF][START_REF] Maker | NIKE3D a non-linear, implicit, three-dimensional finite element code for solid and structural mechanics[END_REF]). In order to achieve good behavior, we feel it is imperative to use some form of "projection" to reduce the number of volumetric constraints. This is absolutely essential for lower-order elements, and very important for higher-order elements as well (as we will show herein). Designing projection schemes for higher-order elements, as far as we are aware, has not been pursued previously in the literature. The B scheme (see Hughes [START_REF] Hughes | Equivalence of finite elements for nearly incompressible elasticity[END_REF]) is a formalism that utilizes projection and we propose a new family of higher-order B schemes in the sequel. Engineering and mathematical points of view are both relevant to this pursuit, but are not entirely in accord. In these cases, we favor the view gained from engineering experience (de gustibus non est disputandum). In particular, some elements produced by our scheme are definitely known not to satisfy the well-known Ladyzenskaya-Babuska-Brezzi (LBB) stability condition (see, e.g., [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Brezzi | On some mixed finite element methods for incompressible and nearly incompressible finite elasticity[END_REF][START_REF] Szabó | Finite Element Analysis[END_REF]), but nevertheless are favored and widely utilized in engineering applications, the prime example being the mean dilatation, bilinear quadrilateral element "Q 1 /Q 0 ". It should be mentioned here that this element fails to satisfy the LBB condition only "slightly" so is not altogether mathematically abhorrent. We also wish to emphasize we are certainly not opposed to satisfying the LBB condition, and some of our other constructs may, at least under certain situations, but we also want to develop methods that have the highest probability of being used in engineering applications. This is the key issue in the approaches we have proposed. The large-deformation counterpart involves projection of the deformation gradient, a so-called F-scheme, involving a product decomposition into volumetric and deviatoric factors. There has been only limited use of such schemes in the literature, and in our opinion none is completely satisfactory. More about this later. A second contribution of this work is the development of a new F formulation based on a modified minimum potential energy principle. We derive the variational equations and a consistent tangent operator for general hyperelastic materials. There are a number of "unexpected" terms in the consistent tangent and it is symmetric, a theoretical improvement over some previous F formulations. We also perform a number of numerical calculations on quasi-static, small-and large-deformation test problems. The calculations are very encouraging and seem to validate the new methods.

Literature survey

Several techniques are available in the literature to deal with incompressibility. The mixed Galerkin finite element formulation is one of the most common. However, the order of approximation for both pressure and displacement/velocity variables cannot be chosen arbitrarily. The method must fulfill the LBB condition to ensure stability and optimal convergence (see, e.g., [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Brezzi | On some mixed finite element methods for incompressible and nearly incompressible finite elasticity[END_REF][START_REF] Szabó | Finite Element Analysis[END_REF]) or must be used in conjunction with stabilization techniques which allow the use of a wide variety of interpolations, including equal-order. These stabilization techniques have been utilized with success in the past twenty years in fluid mechanics; the most commonly used are Streamline Upwind Petrov-Galerkin (SUPG, see, e.g., [START_REF] Brooks | Streamline-Upwind-Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF][START_REF] Hughes | A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuska-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accomodating equal order interpolations[END_REF]) and Galerkin Least Squares (GLS, see, e.g., [START_REF] Hughes | A new finite element formulation for computational fluid dynamics: VIII. The Galerkin Least Squares method for advective-diffusive equations[END_REF]). In more recent work, these ideas have been successfully extended to nonlinear solid mechanics [START_REF] Klass | A stabilized mixed finite element method for finite elasticity. Formulation for linear displacement and pressure interpolation[END_REF][START_REF] Maniatty | Higher order stabilized finite element method for hyperelastic finite deformation[END_REF][START_REF] Ramesh | Stabilized finite element formulation for elastic-plastic finite deformations[END_REF]. Another possible form of stabilization for displacement based formulations is "hourglass" control used in conjunction with reduced integration rules (see, for example, [START_REF] Belytschko | Hourglass control in linear and nonlinear problems[END_REF][START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF][START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF][START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF]).

A different approach to handle the incompressibility condition is the incompatible mode technique proposed by Wilson et al. [START_REF] Wilson | Incompatible displacement models[END_REF], and later generalized as the enhanced or assumed strain method (EAS) by Simo et al. [START_REF] Simo | Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes[END_REF][START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF] within a three-field variational principle. Although this method has exhibited hourglass instabilities in some cases, several improvements have been proposed to overcome these effects, in particular, by using an alternative enhancement strategy [START_REF] Armero | On the formulation of enhanced strain finite elements in finite deformation[END_REF], or a mixed-enhanced strain formulation [START_REF] Kasper | A mixed-enhanced strain method. Part I: geometrically linear problems[END_REF][START_REF] Kasper | A mixed-enhanced strain method. Part II: geometrically nonlinear problems[END_REF].

The approach of greatest engineering interest is to use a simple puredisplacement formulation. This idea goes back in the 1970's with the so-called "reduced" and "selective" integration methods, and B technique [START_REF] Hughes | Generalization of selective integration procedure to anisotropic and nonlinear media[END_REF], for which equivalence theorems with mixed methods have been obtained [START_REF] Hughes | Equivalence of finite elements for nearly incompressible elasticity[END_REF][START_REF] Hughes | A general penalty/mixed equivalence theorem for anisotropic, incompressible finite elements[END_REF][START_REF] Malkus | Mixed finite element methods -Reduced and selective integration techniques: a unification of concepts[END_REF]. If seen in a more general way, as a strain projection technique [START_REF] Hughes | Generalization of selective integration procedure to anisotropic and nonlinear media[END_REF][START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], one can include within the same B framework the well known mean-dilatational formulation of Nagtegaal et al. [START_REF] Nagtegaal | On numerically accurate finite element solutions in the fully plastic range[END_REF]. This formulation is widely used in large-scale and commercial codes [START_REF]LS-DYNA Keyword User's Manual[END_REF][START_REF] Maker | NIKE3D a non-linear, implicit, three-dimensional finite element code for solid and structural mechanics[END_REF]. The major numerical developments on F-type methods, the generalization of B-type methods to hyperelastic finite-deformation formulation, have been done by Hughes et al. [START_REF] Hughes | Finite element formulation and solution of contact-impact problems in continuum mechanics-III[END_REF], Simo et al. [START_REF] Simo | Variational and projection methods for the volume contraint in finite deformation elasto-plasticity[END_REF] (within a three-field Hu-Washizu variational principle), and de Souza Neto et al. [START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] De Souza Neto | F-bar based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF], who proposed an F technique for linear elements. These works exhibit some shortcomings: a consistent tangent was not derived in [START_REF] Hughes | Finite element formulation and solution of contact-impact problems in continuum mechanics-III[END_REF], it was never published in the open litterature, and not pursued further; [START_REF] Simo | Variational and projection methods for the volume contraint in finite deformation elasto-plasticity[END_REF] is a mixed method requiring pressure unknowns; and [START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] De Souza Neto | F-bar based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF] leads to a non-symmetric consistent tangent.

The development of isogeometric analysis, recently proposed by Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], introduced a new numerical method in an attempt to improve ge-ometry, solution representation and mesh refinement compared with finite element analysis. Isogeometric analysis is based on the geometric primitives of computer graphics and Computer-Aided Design (CAD). The first implementation of isogeometric analysis was based on Non-Uniform Rational B-Splines (NURBS), a standard technology in CAD [START_REF] Cohen | Geometric Modeling with Splines: An Introduction[END_REF][START_REF] Farin | NURBS Curves and Surfaces: from Projective Geometry to Practical Use[END_REF][START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF][START_REF] Rogers | An Introduction to NURBS With Historical Perspective[END_REF]. It was first applied to three dimensional linear elasticity and advection-diffusion problems in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]. The ability of NURBS to represent precise geometry was successfully applied to solve problems of vascular fluid-structure interaction with patient-specific geometries in [START_REF] Bazilevs | Isogeometric fluidstructure interaction analysis with applications to arterial blood flow[END_REF][START_REF] Zhang | Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow[END_REF]. The mathematical theory of NURBS based isogeometric analysis was begun in [START_REF] Bazilevs | Isogeometric Analysis: approximation, stability and error estimates for h-refined meshes[END_REF], where stability proofs and optimal error estimates were obtained assuming h-refinement. An interesting feature of isogeometric analysis, not shared by finite element analysis, is so-called k -refinement, which simultaneously increases polynomial order and continuity of the basis functions. Using this feature, isogeometric analysis attains better accuracy than finite elements for structural vibrations [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF]. The higher continuity, beyond the geometric interest, has also proved advantageous for turbulent flows problems [START_REF] Akkerman | The role of continuity in residual-based variational multiscale modeling of turbulence[END_REF]. The above developments and observations suggest NURBS-based isogeometric analysis may offer some new possibilities in large scale structural mechanics applications.

Outline

The paper is organized as follows: In Section 2 the basics of isogeometric analysis are recalled with an emphasis on k -refinement. In Section 3, the use of high-order NURBS within a projection technique is studied in the geometrically linear case with a B method to investigate the choice of approximation and projection spaces with NURBS. The performance of the proposed strategy is demonstrated on well-known numerical examples. Section 4 presents the new F formulation. The consistent tangent operator of the Newton-Raphson iterative scheme for the fully nonlinear elastic case is derived in detail in Appendix B. The performance of the new method is compared with some of the techniques cited earlier on several numerical examples in Section 5, and the good behavior of the new methodology is again observed. Finally, Section 6 presents some conclusions and possible directions of future work.

Brief overview of Isogeometric Analysis

This section gives a brief overview of NURBS-based isogeometric analysis. A more detailed description can be found in [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF][START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]. NURBS functions are a standard tool in CAD to describe and model curves and surfaces [START_REF] Cohen | Geometric Modeling with Splines: An Introduction[END_REF][START_REF] Farin | NURBS Curves and Surfaces: from Projective Geometry to Practical Use[END_REF][START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF][START_REF] Rogers | An Introduction to NURBS With Historical Perspective[END_REF].

B-Splines

NURBS are built from B-Splines. B-Splines are piecewise polynomial curves composed of a linear combination of B-Spline basis functions. The B-Spline parametric space is local to "patches" in contrast with finite elements in which each element carries its own parametrization. Patches may be thought of as subdomains.

Knot Vectors

A knot vector in one dimension is a set of coordinates in the parametric space written Ξ = {ξ 1 , ..., ξ n+p+1 },

where ξ i ∈ R is the i th knot, i is the knot index, i = 1, 2, ..., n + p + 1, p is the polynomial order of the B-Spline and n is the number of basis functions corresponding to it. The knots partition the parameter space into elements, and the interval [ξ 1 , ξ n+p+1 ] constitutes a patch. A knot vector is said to be uniform if its knots are uniformly spaced and non-uniform otherwise. Knot values may be repeated, that is, more than one knot may take a given value. ], but are not, in general, interpolatory at the interior knots. Open knot vectors enable patches to be assembled in essentially the same way elements are assembled in finite element analysis.

Basis functions

B-Spline basis functions are defined recursively starting with piecewise constants (p = 0); given a knot vector Ξ

N i,0 (ξ) =    1 if ξ i ≤ ξ < ξ i+1 , 0 otherwise. ( 2
)
For p > 1, they are defined by:

N i,p (ξ) = ξ -ξ i ξ i+p -ξi N i,p-1 (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 N i+1,p-1 (ξ). (3)
An example for a uniform knot vector is presented in Figure 1. Note that for p = 0 and p = 1, the basis functions are the same as for standard piecewise constant and linear finite element functions. An example of quadratic basis function for an open non-uniform knot vector is presented in Figure 2. Note that the basis functions are interpolatory at the ends of the interval and also at ξ = 4, the location of a repeated knot, where only C 0 -continuity is attained. Basis functions of order p have p -m i continuous derivatives across knot ξ i , where m i is the multiplicity of the value of ξ i in the knot vector. When the multiplicity of a knot is exactly p, the basis is interpolatory at that knot. When the multiplicity is p + 1, the basis is discontinuous and the patch is split into two separate patches.

1 0 0 1 2 3 4 5 N 7,2 N 2,2 N 1,2 N 3,2 N 4,2 N 5,2 N 6,2 N 8,2
ξ Fig. 2. Quadratic basis functions for an open, non-uniform knot vector ξ = {0, 0, 0, 1, 2, 3, 4, 4, 5, 5, 5}

An important property of B-Spline basis functions formed from an open knot vector is that they constitute a partition of unity, that is, ∀ξ:

n i=1 N i,p (ξ) = 1. ( 4
)
This feature is also shared with finite elements and meshless methods. Another interesting property is that the support of each N i,p is compact and contained in the interval [ξ i , ξ i+p+1 ]. Finally, one can note that each basis function is point-wise non-negative over the entire domain of definition: N i,p ≥ 0 ∀ξ. This means that all of the entries of the mass matrix will be positive, which has implications for developing mass lumping schemes (see [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF] for an initiatory investigation).

B-Spline curves

B-Spline curves in R d are constructed by taking linear combinations of B-Spline basis functions. The vector-valued coefficients of the basis functions are referred to as control points. They are analogous to nodal coordinates in finite element analysis in that they are the coefficients of the basis functions. However, the non-interpolatory nature of the basis does not lead to the usual finite element geometric interpretation of the control points. The piecewise linear interpolation of the control points defines the control net. Given n basis functions, N i,p , i = 1, 2, ..., n, and n corresponding control points B i ∈ R d , i = 1, 2, ..., n a piecewise polynomial B-Spline curve is given by:

C(ξ) = n i=1 N i,p (ξ)B i . ( 5
)
An example is shown in Figure 3 for the quadratic basis illustrated in Figure 2. Note that the curve is interpolatory at the first and last control points, a general feature of a curve built from an open knot vector. The curve is also interpolatory at the sixth control point. This is due to the fact that the multiplicity of the knot ξ = 4 is equal to the polynomial order there, that is, m 5 = p = 2. Note also that the curve is tangent to the control polygone at the first, sixth and last control points. The curve is C p-1 = C 1 -continuous everywhere except at the repeated knot ξ = 4 where it is C p-m 5 = C 0 -continuous. The properties of B-Spline curves follow directly from the properties of their basis functions. For example, B-Spline curves have continuous derivatives through order p -1 in the absence of repeated knots or control points. Repeating a control point or a knot k times decreases the number of continuous derivatives by k. An affine transformation of a B-Spline curve is obtained by applying the transformation to the control points. This turns out to be the essential property to guarantee satisfaction of so-called "patch-tests", as discussed in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]. This property is referred to as affine covariance.

2.1.4 h-and p-refinement: knot insertion and order elevation.

The analogue of h-refinement is knot insertion. Knots can be inserted without changing the curve parametrically or geometrically. Given a knot vector Ξ = {ξ 1 , .

B i = α i B i + (1 -α i )B i-1 , (6)
where

α i =          1 if 1 ≤ i ≤ k -p, ξ-ξ i ξ i+p -ξ i if k -p + 1 ≤ i ≤ k, 0 if k + 1 ≤ i ≤ n + p + 2. (7)
Knot values already present in the knot vector can be repeated. However, this reduces the continuity of the basis at the corresponding knot. The continuity of the curve is preserved by choosing the control points using Eqs. ( 6) and [START_REF] Bazilevs | Isogeometric Analysis: approximation, stability and error estimates for h-refined meshes[END_REF]. The mechanism for implementing p-refinement is order elevation. The polynomial order of basis functions can be increased without changing the curve parametrically or geometrically. Note that each unique knot value in Ξ must be repeated in order to preserve discontinuities in the p th derivatives of the curve that is being elevated. The number of new control points depends on the multiplicities of existing knots. As for knot insertion, the solution space spanned by the elevated basis contains the space spanned by the original basis. Thus, it is possible to order elevate without changing the geometry and the parametrization of the B-Spline curve (see [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF][START_REF] Rogers | An Introduction to NURBS With Historical Perspective[END_REF] for further details and algorithms).

An example of h-and p-refinement is presented on Figure 4. The original curve consists of quadratic B-Splines with knot vector Ξ = {0, 0, 0, 1, 1, 1}. The original curve, knot vector and basis functions are shown on the left. A new knot is inserted at ξ = 0.5. The new curve, knot vector and basis functions are shown in the center of the figure. Note that the curve remains unchanged but that basis functions and control points are different; there is one more of each. Next, the curve is order elevated once; the new curve, control points and basis functions can be seen on the right of the figure. This time the multiplicity of the knots is increased by one. The location and number of control points change, but the curve remains the same. There are now four cubic basis functions. Despite the same number of basis functions, note that the locations of the new control points obtained from knot insertion and order elevation are different.

2.1.5 k-refinement: higher-order and higher continuity

The last paragraph presented the analogues of finite element h-and p-refinement: knot insertion and order elevation. To be identical to hrefinement, knot insertion must be performed such that each new knot has a multiplicity equal to the polynomial order of the basis, ensuring C 0 continuity at each knot. Similarly, if we begin with a mesh where all functions are already C 0 , order elevation coincides exactly with the usual notion of p-refinement in the finite element literature. However, knot insertion and order elevation provide us with additional possibilities (see [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF]).

A higher-order alternative elevation strategy present itself; it takes advantage of the fact that knot insertion and order elevation do not commute. If a unique knot value ξ is inserted between two distinct knot values in a curve of order p, the number of continuous derivatives of the basis functions at ξ is p -1.

If we subsequently elevate to order q, the multiplicity of every knot value is increased to maintain C p-1 continuity of the basis functions at ξ. If instead we order elevate the original curve to q and then insert a unique knot value ξ, the basis would have q -1 continuous derivatives at ξ. This strategy is referred to as k-refinement; it has no analogue in standard finite element analysis.

Ξ = {0, 0, 0, 1, 1, 1} p = 2 p = 1 → p = 2
Fig. 5. Non commutativity of h-and p-refinement. The starting point is one linear element. Classical p-refinement strategy: knot insertion followed by order elevation results in seven C 0 quadratic functions. New higher-order k-refinement strategy: order elevation followed by knot insertion results in five C 1 quadratic functions.

The concept of k-refinement is important because isogeometric analysis is fundamentally a higher-order approach. In traditional p-refinement there is a very inhomogeneous structure to arrays due to the different basis functions associated with surface, edge, vertex, and interior nodes. Furthermore, maintaining C 0 continuity during the refinement process implies a proliferation in the number of nodes, In k-refinement, there is a homogeneous structure within patches and growth in the number of control variables is limited (see [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF]).

An example of the k-refinement process and its comparison with traditional prefinement is given in Figure 5. This example shows that k-refinement results in fewer functions with higher continuity, thus fewer control variables or degrees of freedom. Starting with p + 1 functions, inserting n -(p + 1) knots (i.e., to obtain n basis functions), followed by r order elevations, results in (r+1)n-rp C p-1 basis functions. Using k-refinement, that is, starting with p + 1 functions, performing r order elevations, followed by the insertion of n -(p + 1) knots, results in n + r C r+p-1 basis functions. More details on this can be found in [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF][START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF][START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]. Keeping in mind that the preceding numbers are raised to the d power in d dimensions, it is clear that the k-refinement strategy produces fewer unknowns than p-refinement, for the same mesh and order of approximation.

In the k-refinement process on a fixed mesh, knots are added at the boundary points, increasing their multiplicity, but no interior knots are added. It is interesting to note that in the periodic case, in which there are no boundaries and consequently no open knot vectors (see, e.g. Bazilevs et al. [START_REF] Bazilevs | Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[END_REF]), no additional equations are engendered by k-refinement.

B-Spline surfaces and solids

Given a control net {B i,j }, i = 1, ..., n, j = 1, ..., m and knot vectors Ξ = {ξ 1 , ..., ξ n+p+1 }, H = {η 1 , ..., η m+q+1 }, a tensor product B-Spline surface is defined by:

S(ξ, η) = n i=1 m j=1 N i,p (ξ)M j,q (η)B i,j , (8)
where N i,p and M j,q are B-Spline basis functions of order p and q respectively. For purposes of numerical integration of arrays constructed from B-Splines, "elements" are taken to be knot spans, namely [ξ i , ξ i+1 ] × [η j , η j+1 ], as seen previously in the one-dimensionnal case in Figure 3. As presented in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], integrals are pulled back to the parent element by the classical change of variables formula and standard Gaussian quadrature rules are employed as in finite elements (see for example [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], Chapter 3). However, we note that optimally efficient quadrature rules for B-Spline patches do not seem to be known yet.

Analogously, tensor product B-Spline solids can be defined; given a control net {B i,j,k }, i = 1, ..., n, j = 1, ..., m, k = 1, ..., l and knot vectors Ξ = {ξ 1 , ..., ξ n+p+1 }, H = {η 1 , ..., η m+q+1 } and Z = {ζ 1 , ..., ζ l+r+1 }:

S(ξ, η, ζ) = n i=1 m j=1 l k=1 N i,p (ξ)M j,q (η)L k,r (ζ)B i,j,k , (9)
where N i,p , M j,q and L k,r are B-Spline basis functions of order p,q and r respectively.

Rational B-Splines

As described previously, NURBS are built from B-Splines. Specifically, NURBS entities in R d can be built from a projective transformation of B-Spline entities in R d+1 . In particular, conic sections, such as circles and ellipses, can be exactly constructed by projective transformations of piecewise rational quadratic curves.

To obtain a NURBS curve in R d , we start from B w i , i = 1, ..., n, a set of control points ("projective points") for a B-Spline curve in R d+1 with knot vector Ξ. The control points for the NURBS curve are defined by:

(B i ) j = (B w i ) j w i , j = 1, ..., d, (10)
where (B i ) j is the j th component of the vector B i and w i = (B w i ) d+1 is referred to as the i th weight. The rational basis functions and NURBS curve are given by:

R p i (ξ) = N i,p (ξ)w i n i=1 N i,p (ξ)w i , (11)
C(ξ) = n i=1 R p i (ξ)B i . (12)
Rational surfaces and solids are defined analogously in terms of the rational basis functions; see [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF][START_REF] Rogers | An Introduction to NURBS With Historical Perspective[END_REF] for further details:

R p,q i,j (ξ, η) = N i,p (ξ)M j,q (η)w i,j n i=1 m j=1 N i,p (ξ)M j,q (η)w i,j , (13)
R p,q,r i,j,k (ξ, η, ζ) = N i,p (ξ)M j,q (η)L k,r (ζ)w i,j,k n i=1 m j=1 l k=1 N i,p (ξ)M j,q (η)L k,r (ζ)w i,j,k . (14)
In the following, we summarize noteworthy properties of NURBS:

• NURBS basis functions formed from an open knot vector constitute a partition of unity. • The continuity and support of NURBS basis functions are the same as for B-Splines. • NURBS possess the property of affine covariance.

• If all weights are equal, NURBS become B-Splines.

• NURBS surfaces and solids are the projective transformations of tensor product piecewise polynomial entities.

NURBS as a basis for Isogeometric Analysis

The isogeometric analysis framework, based on NURBS consists of the following features: • A mesh for a NURBS patch is defined by the product of knot vectors. For example in three dimensions, a mesh is given by Ξ × H × Z. • Knot spans subdivide the domain into "elements."

• The support of each basis function consists of a small number of elements.

• The control points associated with the basis functions define the geometry.

• The isoparametric concept is invoked. The coefficients of the basis functions are the degrees of freedom or control variables. • Three different mesh refinement strategies are possible: analogues of h-and p-refinement and the new higher-order k-refinement scheme. • Element arrays constructed from isoparametric NURBS are assembled into global arrays in the same way as finite elements (see [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], Chapter 2). • Dirichlet boundary conditions are applied to the control variables. Homogeneous conditions are satisfied pointwise. For inhomogeneous conditions, the boundary values must be approximated by functions lying in the NURBS space; this results in "strong" but approximate satisfaction. Another option is to impose Dirichlet conditions "weakly" (see, for example, [START_REF] Bazilevs | Weak Dirichlet boundary conditions for wall-bounded turbulent flows[END_REF]). Neumann boundary conditions are satisfied naturally. The situation is very similar to standard finite elements (see [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], Chapters 1 and 2).

Table 1 summarizes the similar and different properties of isogeometric analysis and finite element analysis.

3 B formulation for linear elasticity using NURBS

The boundary value problem of compressible and incompressible linear elasticity

The boundary value problem of compressible elasticity for a body Ω is given by the following:

Given f : Ω → R 3 , g : Γ g → R 3 , and h : Γ h → R 3 , find u : Ω → R 3 such that: divσ + f = 0 in Ω, ( 15
) u = g on Γ g , ( 16
) σ • n = h on Γ h , (17)
n is the exterior unit normal on Γ, the boundary of Ω, g is the prescribed displacement on Γ g and h is the prescribed traction on Γ h , which form together the boundary Γ = Γ h ∪ Γ g of Ω, and f is the body force. The stress tensor σ is defined in terms of the strain tensor ε by the generalized Hooke's law:

ε = ∇ s u = 1 2 (∇u + ∇u T ) or ε ij = 1 2 ∂u i ∂x j + ∂u j ∂x i , (18)
σ = c : ε or σ ij = c ijkl ε kl (19)
The Einstein summation convention is employed for spatial indices throughout.

In the compressible isotropic linear elastic case, Hooke's law can be expressed in terms of the Lamé parameters λ and µ by:

c ijkl = λδ ij δ kl + µ(δ ik δ jl + δ il δ jk ), (20)
σ ij = λu k,k δ ij + 2µε ij , (21)
where

λ = ν2µ (1 -2ν) , µ = E 2(1 + ν) , ( 22
)
and ν is Poisson's ratio and E is Young's modulus.

Clearly, as ν → 1 2 , λ approaches infinity. The value ν = 1 2 thus represents incompressibility. The constitutive equation needs to be modified in this case:

σ ij = -pδ ij + 2µε ij , (23)
where p, the hydrostatic pressure, is determined as part of the solution of the boundary value problem. As p represents an additional unknown, the kinematic condition of incompressibility must be introduced as an additional equation:

divu = u k,k = 0 in Ω. ( 24
)
However, to represent the nearly incompressible case, that is, the ratio λ/µ is large but not infinite, the compressible theory can be applied with some modifications in the discrete case.

Selective and reduced integration

Among the original techniques to successfully solve the nearly incompressible case are selective and reduced integration techniques as presented in [START_REF] Hughes | Equivalence of finite elements for nearly incompressible elasticity[END_REF][START_REF] Hughes | Generalization of selective integration procedure to anisotropic and nonlinear media[END_REF][START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF].

The basic idea of selective integration is to split the bilinear form of the variational formulation into its λ and µ contributions. Using vector representation of tensors, and introducing the strain-displacement matrix B (see [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF], chapter 2 for a definition of B), the element stiffness matrix is:

k e = Ω e B T DBdΩ. (25)
The material properties matrix D can be split into its λ and µ parts:

D = µ D + λ D, (26)
which results in the same split for the element stiffness matrix:

k e = µ ke + λ ke . (27)
Since λ/µ ≫ 1, the numerical values of the entries in the second matrix in Eq. ( 27) tend to be very large compared with those of the first. The idea is to use a standard Gauss quadrature rule to compute ke and a reduced (i.e., lower order Gauss quadrature) rule on ke to reduce the number of incompressibility constraints. The simplest example is the piecewise bilinear quadrilateral element in plane strain for which the normal rule is the 2 × 2 Gauss points rule and the reduced rule is the one-point Gauss rule. Equivalence theorems between reduced/selective integration elements and mixed formulation elements have been obtained by Hughes and Malkus [START_REF] Hughes | Equivalence of finite elements for nearly incompressible elasticity[END_REF][START_REF] Hughes | A general penalty/mixed equivalence theorem for anisotropic, incompressible finite elements[END_REF][START_REF] Malkus | Mixed finite element methods -Reduced and selective integration techniques: a unification of concepts[END_REF]. These theorems show that selective and reduced integration techniques are a simple way of attaining the performance of mixed methods without engendering unwanted degrees of freedom. This seems to be a primary reason why these techniques are used in large-scale and commercial codes [START_REF]LS-DYNA Keyword User's Manual[END_REF][START_REF] Maker | NIKE3D a non-linear, implicit, three-dimensional finite element code for solid and structural mechanics[END_REF].

Strain projection and the B method

Some drawbacks of the selective and reduced integration rules are that the equivalence theorems with mixed formulations are not valid in the axisymmetric case and that these rules are difficult to generalize in anisotropic cases. This can become a particular problem with nonlinear problems since the tangent moduli always exhibits anisotropic character. These difficulties were overcome by introducing a strain projection approach, referred to as the B method (Hughes [START_REF] Hughes | Generalization of selective integration procedure to anisotropic and nonlinear media[END_REF]).

The main idea in the strain projection approach is to additively split the strain tensor into its deviatoric and dilatational (i.e., volumetric) parts

ε(u) = ε dev (u) + ε dil (u), (28)
where

ε dil (u) = 1 3 (divu) I or ε dil ij (u) = 1 3 ∂u k ∂x k δ ij , ( 29
)
and I is the identity tensor.

To achieve an effective formulation in the nearly incompressible case, the dilatational part is replaced by an "improved" dilatational contribution (e.g., a "projected" one), using a linear projection operator π

εdil (u) = π ε dil (u) . (30)
In terms of B, the strain-displacement matrix, we have

B = B dev + Bdil , (31)
with

B dev = B -B dil . ( 32
)

Weak form of the B method

Let us recall the usual principle of minimum potential energy to obtain the weak form in linear elasticity. We start by defining the trial and weighting spaces

S = {u | u ∈ H 1 (Ω), u| Γg = g} and V = {w | w ∈ H 1 (Ω), w| Γg = 0}.
Given u ∈ S the potential energy is defined by:

Π(u) = Ω Ψ(ε(u))dΩ - Ω u • f dΩ - Γ h u • hdΓ, ( 33
)
where Ψ is the strain energy density:

Ψ(ε(u)) = 1 2 c ijkl ε ij (u)ε kl (u). ( 34
)
Minimizing the potential energy is equivalent to satisfying the variational equation, or weak form; that is given u ∈ S, w ∈ V and ǫ ∈ R,

∂Π(u + ǫw) ∂ǫ | ǫ=0 = 0 ⇔ a(w, u) = (w, f ) + (w, h) Γ h , (35)
where

a(w, u) = Ω ε ij (w)c ijkl ε kl (u)dΩ, ( 36
) (w, f ) = Ω u • f dΩ, (37)
(w, h) Γ h = Γ h u • hdΓ. ( 38
)
Considering now the potential energy in terms of the modified strain tensor ε:

Π(u) = Ω Ψ(ε(u))dΩ - Ω u • f dΩ - Γ h u • hdΓ, (39)
and recalling that ε is a linear operator, one can write:

ε(u + ǫw)) = ε(u) + ǫε(w). (40)
Therefore, the directional derivative of the modified strain energy density is given by:

∂Ψ(ε(u + ǫw)) ∂ǫ | ǫ=0 = ∂Ψ ∂ε ij (ε(u)) εij (w), (41)
which allows us to define the modified stress tensor as:

σij = c ijkl εkl . ( 42
)
The minimization of the modified potential energy finally becomes:

0 = Ω εij (w)c ijkl εkl (u)dΩ - Ω w • f dΩ - Γ h w • hdΓ, (43)
which represents the variational equation for the B method. The first integral in Eq. ( 43) is a bilinear form, and expressing it the following way,

ā(w, u) = Ω εij (w)c ijkl εkl (u)dΩ, (44)
the B variational formulation can be stated as:

Find u ∈ S, such that ∀w ∈ V ā(w, u) = (w, f ) + (w, h) Γ h . (45)
3.4 B and higher-order NURBS

Projection operator and space with NURBS

The use of the B method within isogeometric analysis requires further investigation into the choices of the projection operator and the associated space onto which the projection will be performed. Since the technique has been applied mostly to piecewise bilinear and trilinear finite elements, and we want to make intensive use of the properties of high-order k -refined NURBS, these topics need to be studied without any assumption on the order of approximation.

In the discrete case, we have

u h (x) = n A=1 u A N A (x), (46)
likewise

w h (x) = n A=1 w A N A (x), (47)
where N A are the NURBS basis functions and u A and w A are the associated control variables. This reduce to the usual shape functions and degrees of freedom when classical finite elements are utilized.

In developing the B method for higher-order finite elements and NURBS, we need to define the linear projection operator and the spaces upon which to project the dilatational strain. Throughout, we use the L 2 projection of the strains. For the spaces, we define the following procedure: assume the displacement space is given. We shall refer to it as Q p , that is quadrilateral, or hexahedral, elements of order p. The continuity of Q p elements within a patch can be any order k from 0 to p-1. In this paper, we are particularly interested in elements of maximal continuity, namely, C k = C p-1 . We always assume an open knot vector construction so that only C 0 continuity is attained across patch interface. The basis functions for the projected dilatational strain are taken to be one order lower, and usually one order of continuity lower, namely the space Q p-1 , of continuity C k = C p-2 . The only exception occurs when p ≥ 2, but there are lines or surfaces of C 0 continuity within a patch. In this situation, the projected space is also taken to be C 0 continuous across those lines or surfaces. There is nothing fundamental about this choice. It is simply a convenience due to the data structure we employ in our code. In practice, C 0 lines or surfaces within a patch exist when conic sections such as circles, are constructed by standard NURBS algorithms (see e.g., Piegl and Tiller [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF] and Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]).

As an example of what our construction produces, consider the displacement space Q 1 . This is the space of bilinear quadrilaterals, or trilinear hexahedra, and is C 0 continuous across element boundaries (which correspond to knots in this case). The space for projected dilatational strain is then Q 0 , of continuity class C -1 , that is, piecewise constants. This element becomes the classical mean dilatational element (see Hughes and Allik [START_REF] Hughes | Finite elements for compressible and incompressible continua[END_REF], Nagtegaal et al. [START_REF] Nagtegaal | On numerically accurate finite element solutions in the fully plastic range[END_REF] and Hughes [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]) referred to, herein, as Q 1 /Q 0 . We wish to emphasize that we recognize that there is no guarantee that our construction produces discrete approximations satisfying the LBB condition (see e.g., [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Brezzi | On some mixed finite element methods for incompressible and nearly incompressible finite elasticity[END_REF][START_REF] Szabó | Finite Element Analysis[END_REF]). In fact, it is entirely obvious that this is the case, the Q 1 /Q 0 element being a prime example of an element that fails. However, despite the theoretical deficiency, it must be noted that Q 1 /Q 0 is without doubt the most utilized element in the history of nonlinear structural finite element analysis, and it is widely used in fluids as well. There are other practical issues at play here (see Hughes [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]). We will say, however, that our construction produces a good balance between the number of displacement degrees of freedom and number of dilatational constraints on each patch. Asymptotically, "constraint ratios" are 2 in two dimensions and 3 in three dimensions, in all cases (see Hughes [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF] for a discussion of constraint ratios). But from the LBB perspective, we may have slightly too many constraints, at least for some cases. We note that for the discretizations of main interest herein, namely p ≥ 2, Q p /Q p-1 , of internal patch continuity C p-1 /C p-2 , there are no known theoretical results, one way or the other, at this point in time. However, results from the p-method community seem to indicate that incompressible locking is not an issue for sufficiently great p (see Suri [START_REF] Suri | Analytical and computational assessment of locking in the hp finite element method[END_REF], Szabó and Babuska [START_REF] Szabó | Finite Element Analysis[END_REF] and Szabó et al. [START_REF] Szabó | The p-version of the finite element method[END_REF]).

The construction we advocate takes a middle road. It does not reduce the number of constraints sufficiently to satisfy the LBB condition for all cases, but it does reduce it to a value that seems to attain a good balance in all cases. Our experience with it is good, as the numerical results will indicate in the sequel. On the other hand, our experience with higher-order elements without projection is that they are overconstrained and not competitive with the reduced constraint, projection elements advocated. Likewise, spectral elements of order Q p /Q p-2 , and continuity class C 0 /C -1 are known to satisfy the LBB condition [START_REF] Karniadakis | A spectral element method for solution of two-and three-dimensional time dependent Navier-Stokes equations[END_REF], but convergence is suboptimal and the minimum-order element is Q 2 /Q 0 . In fairness, the spectral element approach is focused on very high-order elements and high precision applications exclusively, not the low-order end of the spectrum. Nevertheless, in nonlinear structural analysis, low-order elements are extremely important. We will now give the technical details of the construction of the spaces. On each patch, open knot vectors are employed. Tensor product constructs are utilized so we focus on the situation in each direction separately. The B-Spline basis is completely determined from Eqs. ( 2) and (3). We need to specify the order of the space and the knot vector. We begin with the displacement space, assumed to be of order p ≥ 1. The knot vector, denoted Ξ p , is assumed to have the following form,

Ξ p = {0, 0, ..., 0 p+1 copies , Ξ int , 1, 1, ..., 1 p+1 copies }, (48)
where, for simplicity, we have assumed the initial and final knots are located at 0 and 1, respectively. The multiplicities of the initial and final knots are p + 1. Ξ int denotes the vector of internal knots. Each internal knot may have a different multiplicity, allowable values being 1, 2, 3, ..., p. The case we are primarily concerned with in this paper is each internal knot having multiplicity 1 which results in maximal smoothness of continuity class C p-1 on each patch. The corresponding knot vector for the projected space, denoted Ξ p-1 , is given by Ξ p-1 = {0, 0, ..., 0

p copies , Ξ int , 1, 1, ..., 1 p copies }. ( 49
)
The order of the projected space is taken to be p -1 ≥ 0. Note that the internal knot vector is exactly the same as for the displacement space, whereas the initial and final knots have multiplicity p. The span of the projected space is precisely the span of the derivatives of all functions in the displacement space. An example of the spaces Q p /Q p-1 , p = 1, 2, 3, in the general case, for 21 a one-dimensional patch of four elements, is given in Figure 6. We see that C p-1 /C p-2 continuity is achieved in all cases. The exception to the general case occurs when there are lines or surfaces of only C 0 continuity within a patch, as mentioned previously. Let us assume that Ξ int has one or more knots having multiplicity p, signifying C 0 continuity. Then, in the space Ξ p-1 , Ξ int needs to be replaced with Ξint , which is identical to Ξ int except for the knots having multiplicity p; in Ξint these knots have multiplicity p -1, preserving C 0 continuity of the projected space within each patch. An example of the spaces Q p /Q p-1 , p = 2, 3, with a point of C 0 continuity inside the patch, is given in Figure 7. Only C 0 continuity is achieved across the repeated knot ξ = 0.75, and C p-1 /C p-2 continuity is achieved elsewhere on the patch interior.

The construction of higher-dimensional B-Splines proceeds by assigning the polynomial order and knot vector for each dimension. For example, in three dimension, we have, Ξ p × H q × Z r and Ξ p-1 × H q-1 × Z r-1 , where

H q = {0, 0, ..., 0 q+1 copies , H int , 1, 1, ..., 1 q+1 copies }, (50)
H q-1 = {0, 0, ..., 0 q copies , H int , 1, 1, ..., 1 q copies }, (51)
Z r = {0, 0, ..., 0 r+1 copies , Z int , 1, 1, ..., 1 r+1 copies }, (52)
Z r-1 = {0, 0, ..., 0 r copies , Z int , 1, 1, ..., 1 r copies }, (53)
and the internal knots in H int and Z int may take multiplicities in the range 1, 2, 3, ..., q, and 1, 2, 3, ..., r, respectively. As before, if there are knots in H int and Z int having multiplicity q and r, then H int and Z int need to be replaced by Hint and Zint , constructed in a similar fashion as to Ξint .

In practice, the N A 's are constructed as

N A (x) = N A • φ -1 (x), (54)
where { N } is a NURBS basis in parametric coordinates and φ is the exact mapping between the parametric and geometric domains. We construct the "tilde basis" { Ñ }, which corresponds to the projection space, by

ÑA (x) = NA • φ -1 (x), (55)
where { N } is the lower order NURBS basis built on the same parametric domain. As described previously, we take the NA 's to be one order lower than the N A 's to reduce the number of incompressibility constraints.

Note that, even in the case of lowest-order elements (i.e., bilinear and trilinear), we still use the exact geometrical mapping. This means our lowest-order elements are isogeometric and precisely fit curved boundaries. We believe that Barth [START_REF] Barth | Numerical methods for gasdynamic systems on unstructured meshs[END_REF] was the first to use this approach, and to demonstrate its effectiveness in compressible fluid calculations.

In the discrete case, Eq. ( 30) becomes:

εdil ij (u h ) = ñ A=1 ÑA εA ij , (56)
where

εA ij = ñ B=1 M -1 AB ÑB , ε dil ij (u h ) Ω = ñ B=1 M -1 AB Ω ÑB ε dil ij (u h )dΩ, (57)
that is

εdil ij (u h ) = ñ A,B=1 n C=1 ÑA M -1 AB Ω ÑB ∂N C ∂x k dΩ u C k δ ij , (58)
and M is the "mass" matrix of the tilde basis, namely

MAB = ÑA , ÑB Ω = Ω ÑA ÑB dΩ. (59)
In summary, the procedure corresponds to L 2 projection of ε dil ij onto the { Ñ } basis.

Numerical consistency

Following the case of the enhanced assumed strain method, one needs to verify that the proposed projection method satisfies the patch test in order to ensure numerical consistency. This property is attained when the average of the projected/assumed strain is the same as the average of the original strain (see, e.g., Simo and Hughes [START_REF] Simo | On the variational foundations of assumed strain methods[END_REF]). Multiplying Eq. ( 58) by the tilde basis and integrating over the domain, one obtains:

Ω ÑD εdil ij (u h )dΩ = ñ A,B=1 Ω ÑD ÑA dΩ M -1 AB Ω ÑB ε dil ij (u h )dΩ, (60)
which can be simplified as:

Ω ÑD εdil ij (u h )dΩ = Ω ÑD ε dil ij (u h )dΩ. ( 61
)
The consistency follows immediately from Eq. ( 61) because the "tilde" lower order NURBS basis is constructed from an open knot vector and therefore constitutes a partition of unity, as presented in Eq. ( 4) in Paragraph 2.1.2. Focusing on this particular point, the main difference between EAS methods and the proposed one is that the tilde functions are not associated to element internal degrees of freedom. In the general case, they are constructed exactly as the NURBS displacement functions (but of one order lower) and therefore share the same properties, such as the partition of unity.

Consequently, any element Q p /Q p-1 , p ≥ 1 is able in particular to represent constant strains and hence the patch test. Moreover, the order of accuracy of the approach is preserved as will be seen in the numerical examples.

Implementational aspects

We consider aspects of solving the global matrix system in this section.

For an isotropic homogeneous linear elastic material, with Hooke's law given by Eq. ( 21), the discrete version of the bilinear form of Eq. ( 44) is given by:

ā(w h , u h ) = n A,B=1 w A i Ω N A ,j ĉijkl N B ,l dΩ u B k + 1 3 (3λ + 2µ) ñ C,D=1 w A i (N A ,i , ÑC ) M -1 CD ( ÑD , N B ,k ) u B k , (62)
where

ĉijkl = µ δ ik δ jl + δ il δ jk - 2 3 δ ij δ kl . (63)
Due to the inverse of M , the second term in Eq. ( 62) increases the population of the stiffness matrix on each patch for p ≥ 2. Note, we always assume use of patches constructed with open knot vectors. This means that the displacement field is continuous across patch interfaces, but no smoother. Consequently, the tilde basis will be discontinuous across patches and M -1 will be uncoupled from patch to patch. Nevertheless, if we use a direct solver to solve the global equation system, we need to account for increased coupling of the equations due to M -1 on each patch. There are at least two ways to circumvent the effect of the increased coupling. One is to use an iterative strategy that does not require the assembly of the stiffness matrix, such as Conjugate Gradients, to solve the global problem. Within each conjugate gradient iteration, a direct solver can be used to evaluate M -1 patch-wise, retaining its sparse band-profile structure. This procedure can be used to solve very large problems. We have used it extensively in our calculations and found it to be very efficient. A second possibility is to replace M with a diagonal, or "lumped" approximation. This would only need to be done in the left-handside matrix, and so would be interpreted as a preconditioner. In this case, the band-profile structure of the preconditioner would be only slightly larger than for the system constructed without projection. Using the consistent M on the right-hand-side would ensure the full accuracy of the projection procedure.

Convergence would require one or more iterations, but this involves only a forward reduction and back substitution for each additional iteration with an existing factorized array when employing a direct solver.

Numerical Examples

We consider NURBS displacement/projection spaces Q p /Q p-1 with p = 1,2,3,4. All internal knots have multiplicity 1, resulting in maximal smoothness, that is C p-1 /C p-2 . For fixed p, we examine "h-refinement," that is, we subdivide the mesh by knot insertion. For a fixed mesh, we also consider k-refinement, that is, we increase p before inserting knots to obtain the desired mesh. k-refinement produces very few additional degrees of freedom (i.e., control variables) as no internal knots are added and only the multiplicities of the initial and final knots are increased by 1. This results in one layer of additional degrees of freedom per dimension. For example, in a one-patch three-dimensional mesh of 3n 3 degrees of freedom, one level of k-refinement results in (3n + 3) 3 degrees of freedom (see Cottrell et al. [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF] for further details).

In the following we present two numerical examples for nearly incompressible linear elasticity to demonstrate the proposed approach. Both are plane strain two-dimensional problems but are solved using a three-dimensional formulation with two linear functions in the out-of-plane direction and boundary conditions enforcing the plane strain constraint. This problem has been solved by many authors to test nearly incompressible formulations under combined bending and shear (see, e.g., [START_REF] Chavan | Locking-free finite element methods for linear and nonlinear elasticity in 2D and 3D[END_REF][START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF][START_REF] Cesar De Sa | New enhanced-strain elements for incompressible problems[END_REF][START_REF] Kasper | A mixed-enhanced strain method. Part I: geometrically linear problems[END_REF]). A tapered panel is clamped on one side and subjected to a uniform shear load on the opposite side. The other sides are traction free. The geometry, loading and boundary conditions are given in Figure 8. The quantity of interest is the vertical displacement of the top right corner of the plate as presented in Figure 8. Its convergence as a function of the number of elements per edge is usually used as a criterion in investigating performance. The results are plotted in Figure 9 for the four orders of approximation with the standard Galerkin approximation, without any treatment of the incompressibility condition, and with the proposed B method. We observe, as expected, that the piecewise bilinear quadrilateral element suffers from severe locking. Increasing the polynomial order of the NURBS functions improves the results, however, we can see that C 1 quadratic NURBS still suffer from locking for coarse meshes. The use of B significantly improves the results, and good convergence is attained even with Q 1 /Q 0 . Note that the results obtained with Q 1 /Q 0 are better than Q 2 and that Q 2 /Q 1 gives better results than Q 4 . Finally, we observe that with B Q 4 /Q 3 we obtain converged results with only a mesh of four elements.

It is also interesting to study convergence as a function of the number of unknowns. The results are shown in Figure 10, where a logarithmic scale has been used for the abscissa. Note here, that in terms of unknowns, Q 1 /Q 0 performs almost the same as Q 3 and Q 4 . We also see that the results with Q 2 /Q 1 and Q 3 /Q 2 are virtually equivalent in terms of unknowns and the best performance is attained with Q 4 /Q 3 . These results demonstrate the superiority of high-order k -refined NURBS: increasing the polynomial order and the continuity only slightly increases the number of unknowns, which enables us to obtain very good results with only a few C 3 quartic functions: for Q 4 /Q 3 , 4 elements and 36 degrees of freedom is enough to obtain the reference value, whereas we need at least 1024 elements and 1089 degrees of freedom to obtain about the same results with Q 1 /Q 0 .

Infinite plate with circular hole under in-plane tension

The next problem is a plane-strain infinite plate with a hole under tension. This problem has been studied previously with isogeometric analysis in [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] assuming an isotropic compressible linear elastic medium. It is interesting from the geometrical point of view because quadratic NURBS can exactly represent the circular hole and the existence of an analytical solution allows us to focus on the convergence rates that the proposed method can attain without geometrical approximation. In the nearly incompressible regime, this problem has been of less interest than the previous one. Some results in the incompressible limit using meshless methods can be found in [START_REF] Dolbow | Volumetric locking in the element free Galerkin method[END_REF][START_REF] Huerta | Locking in the incompressible limit for the Element-Free Galerkin method[END_REF]. The infinite plate is modeled by a quarter plate. The exact solution (Timoshenko and Goodier [START_REF] Timoshenko | Theory of Elasticity[END_REF], pp. 90-93) is evaluated at the boundary of the quarter plate and applied as a Neumann boundary condition. The exact solution is given by the following:

σ rr (r, θ) = T x 2 (1 - R 2 r 2 ) + T x 2 (1 -4 R 2 r 2 + 3 R 4 r 4 ) cos 2θ, ( 64
)
σ θθ (r, θ) = T x 2 (1 + R 2 r 2 ) - T x 2 (1 + 3 R 4 r 4 ) cos 2θ, ( 65
)
σ rθ (r, θ) = - T x 2 (1 + 2 R 2 r 2 -3 R 4 r 4 ) sin 2θ. ( 66
)
The geometry, loading, boundary conditions and parameters are shown in Figure 11. R is the radius of the hole, L the length of the finite quarter plate, and T x is the magnitude of the applied stress for the infinite case. The value for ν was chosen to be very close to the incompressible case to study the convergence of the B isogeometric analysis in that case. A rational quadratic basis is the minimum order capable of representing the geometry. The sequence of meshes obtained from h-refinement used for the convergence study are shown in Figure 12. The geometry of the coarsest mesh can be exactly represented by quadratic NURBS. It consists of two elements, shown on the left of Figure 12, and is defined by the knot vectors Ξ 2 and H 2 , from which the knot vectors Ξ 1 and H 1 are constructed to build the corresponding linear basis:

Ξ 2 × H 2 = {0, 0, 0, 0.5, 1, 1, 1} × {0, 0, 0, 1, 1, 1}, ( 67
) Ξ 1 × H 1 = {0, 0, 0.5, 1, 1} × {0, 0, 1, 1}. ( 68
)
We define the relative error as the error normalized by the corresponding value of the exact solution. Convergence results for the relative error in the L 2 norm of displacement, relative error in the L 2 norm of stress, and relative error in the energy norm are shown in Figure 13. The cubic and quartic NURBS are obtained from k -refinement of the coarsest quadratic mesh. The mesh parameter h is defined as the maximum distance, in the physical space, between diagonally opposite knot locations. As can be seen, the B method obtains good convergence rates in the three norms with relatively coarse meshes. Note that the standard displacement based formulation performs relatively poorly, especially in the L 2 norm of stress, and needs relatively fine meshes to attain convergence rates equivalent to what is obtained with the B formulation. Even when seemingly optimal asymptotic rates of convergence are attained in the case of the standard Q p elements, the error is orders of magnitude greater than for the corresponding projected Q p /Q p-1 B elements. For example, in stress the difference is four orders of magnitude! This result clearly shows that optimal rate of convergence is not the only issue to be considered.

F formulation for nonlinear elasticity

Constitutive equations

The central idea in the development that we propose is, as done in the geometrically linear case, to split the tensor that measures the deformation into its deviatoric (volume preserving) and volumetric-dilatational parts. In the finite deformation case, the deformation gradient F is the relevant tensor, and, contrary to the small deformation case, the split is multiplicative rather than additive. This multiplicative decomposition has been exploited previously by Flory [START_REF] Flory | Thermodynamic relations for highly elastic materials[END_REF], Hughes et al. [START_REF] Hughes | Finite element formulation and solution of contact-impact problems in continuum mechanics-III[END_REF], Simo et al. [START_REF] Simo | Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithm[END_REF][START_REF] Simo | Variational and projection methods for the volume contraint in finite deformation elasto-plasticity[END_REF] (within a three field Hu-Washizu principle), and more recently by de Souza Neto et al. [START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] De Souza Neto | F-bar based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF] in an alternative F approach. The work we present here shares features with these techniques. We want it to be a simple pure displacement formulation, but having greater generality and a more rigorous theoretical base than previous approaches.

We start by writing the multiplicative split of the deformation gradient:

F = F dil F dev , (69)
where detF = J = detF dil and detF dev = 1, (70) which leads to:

F dev = J -1/3 F and F dil = J 1/3 I. ( 71
)
We now define a modified deformation gradient F in terms of the deviatoric part of the deformation gradient F dev and a modified dilatational part of the deformation gradient Fdil :

F = Fdil F dev , (72)
where

Fdil = π(F dil ) = π(J 1/3 )I = J 1/3 I, (73)
with π a linear projection operator identical to the one proposed previously for the linear case. We can finally write, keeping the "bar" notation for the projection:

F = αF, (74)
α = J 1/3 J 1/3 . (75)
We assume hyperelastic homogeneous material behavior for which there exists a free-energy function3 Ψ that depends on the Cauchy-Green tensor C = F T F and from which the second Piola-Kirchhoff stress tensor is derived as:

S = 2 ∂Ψ(C) ∂C (76)
The standard additive decomposition of Ψ (see for example [START_REF] Simo | Computational Inelasticity[END_REF]) into a volumetric part depending only on J and an isochoric part is used:

Ψ(J, C) = Ψ dil (J) + Ψ iso (J, C) (77)

Variational equation and weak form

Let B be the reference configuration of a body and B ′ be its deformed configuration. We introduce a mapping φ : B → B ′ that takes a point X ∈ B to a point x = φ(X) ∈ B ′ . The displacement of a particle from its initial position X to its current position x is given by

u(X) = φ(X) -X = x -X. (78)
The deformation gradient is the second order tensor defined by:

F(X) = ∇ X φ(X) = ∂φ(X) ∂X . (79)
The boundary value problem for finite deformation elasticity for a body with reference configuration B is given by the following:

Given F : B → R 3 , φ : ∂ φ B → R 3 , and H : ∂ H B → R 3 , find φ ∈ S such that: ∇ X • P + F = 0 in B, (80)
φ = φ on ∂ φ B, (81)
P • N = H on ∂ H B, (82)
where N is the unit exterior normal on ∂B, the boundary of B, φ is the prescribed deformation on ∂ φ B and H is the prescribed Piola traction on ∂ H B, which form together the boundary

∂B = ∂ φ B ∪ ∂ H B.
F is the body force per unit initial volume, S is the space of admissible deformations defined by

S = {φ : Ω → R 3 | det(∇ X φ) > 0, φ| ∂ φ B = φ}, ( 83
)
and P is the first Piola-Kirchhoff stress tensor:

P = FS. (84)
Corresponding to what was done previously in the geometrically linear case, given the mapping φ ∈ S, the potential energy is given by

Π(φ) = B Ψ(E(φ)) dV - B φ • F dV - ∂ H B φ • H dΓ, (85)
where E = 1 2 (F T F-I) is the Green-Lagrange strain tensor. The stationarity of the potential energy leads to the variational equation for the finite deformation case (see for example [START_REF] Simo | Computational Inelasticity[END_REF]).

We now introduce the modified potential energy

Π(φ) = B Ψ( Ē(φ)) dV - B φ • F dV - ∂ H B φ • H dΓ, (86)
where the modified Green-Lagrange strain tensor is defined in terms of the modified deformation gradient:

Ē = 1 2 ( FT F -I). (87)
We write the stationarity of the modified potential energy:

∂ Π(φ + ǫW) ∂ǫ | ǫ=0 = 0, (88)
where φ ∈ S,

W ∈ V = {W : B → R 3 | W| ∂ φ B = 0} and ǫ ∈ R.
We need to compute the directional derivative of the strain energy density, keeping in mind that the projection operator is linear:

∂Ψ( Ē(φ + ǫW)) ∂ǫ | ǫ=0 = ∂Ψ(E) ∂E IJ ( Ē(φ)) • ∂ ĒIJ (φ + ǫW) ∂ǫ | ǫ=0 (89)
The calculation of the directional derivative of the modified Green-Lagrange strain tensor is presented in Appendix A, the result is:

∂ Ē(φ + ǫW) ∂ǫ | ǫ=0 = D ĒIJ (φ) • W = 1 3 π(J 1/3 ∇ x • W) J 1/3 -∇ x • W FT F + α FT ∇ X W , ( 90
)
where α is defined in Eq. (75).

We then define the modified second Piola-Kirchhoff stress tensor

SIJ (φ) = S IJ ( Ē(φ)) = ∂Ψ(E) ∂E IJ ( Ē(φ)), (91)
and Eq. (89) becomes

SIJ D ĒIJ (φ) • W = SIJ 1 3 π(J 1/3 ∇ x • W) J 1/3 -∇ x • W FT F + α( FT ∇ X W) IJ . (92)
where

∇ x • W = tr(∇ X WF -1 ) (see Eq. (B.15) in Appendix B).
Defining the modified first Piola-Kirchhoff stress tensor as

PiI = FiJ SIJ , (93)
and introducing the following notation

W i,I = αW i,I + 1 3 π(J 1/3 ∇ x • W) J 1/3 -∇ x • W FiI , (94)
we rewrite Eq. ( 92) as

SIJ D ĒIJ (φ) • W = PiI W i,I . (95)
We introduce the modified Cauchy stress tensor, and push forward to the current configuration Eq. ( 95)

σ = 1 J FS FT , (96)
PiI W i,I = J σij w i,j , (97)
where

w i,j = w i,j + 1 3 π(J 1/3 ∇ x • w) J 1/3 -∇ x • w δ ij (98)
Note that, for all geometrical transformations, we use J = detF, not det F.

The symmetrical part of Eq. ( 98) can be written in direct notation:

∇ x w sym = (∇ x w) sym + 1 3 π(J 1/3 ∇ x .w) J 1/3 -∇ x .w I, = ε(w) + 1 3 (tr [ε(w)] -tr [ε(w)]) I, = ε(w) + εdil (w) -ε dil (w), = ε(w). (99)
Thus, we can write the stationarity of the modified potential energy given in Eq. ( 88) pushed forward to the deformed configuration as

0 = B ′ σij (φ)ε ij (w)dv - B ′ w i f i dv - ∂ h B ′ w i h i dγ, (100)
where f • φ = F, and (h • φ)dγ = HdΓ. This is the finite deformation counterpart of Eq. ( 43), and leads to the corresponding variational formulation.

Find φ ∈ S, such that ∀w ∈ V φ ā(w, φ) = (w, f ) + (w, h) ∂ h B ′ , (101)
where ā(w, φ) is the first integral in Eq. ( 100) and

V φ = {w • φ : Ω → R 3 | w • φ| ∂ φ B = 0}
is the tangent space to S at φ.

Linearized operator

To solve the resultant nonlinear problem, we derive the corresponding consistent linearized operator (see Hughes and Pister [START_REF] Hughes | Consistent linearization in mechanics of solid and structures[END_REF] and Simo and Hughes [START_REF] Simo | On the variational foundations of assumed strain methods[END_REF]) that will be used in the iterative Newton algorithm. Taking inspiration from the linearized operator for the standard compressible finite deformation case and the geometrically linear B theory, we might expect to attain the usual linearized operator in terms of "barred" quantities, namely,

B ′ (w i,j c ijkl u k,l + w k,j σ ij u k,i ) dv (102)
However, as observed by Simo et al. [START_REF] Simo | Variational and projection methods for the volume contraint in finite deformation elasto-plasticity[END_REF], we also find that the structure of the tangent operator associated with the projection method in the fully nonlinear case leads to additional "unexpected" terms.

To obtain the tangent operator we need to compute the directional derivative of the internal virtual work given by the first integral in Eq. (100), or its reference configuration version given in Eq. ( 92):

B ∂ SIJ (φ + ǫU)D ĒIJ (φ + ǫU).W ∂ǫ | ǫ=0 dV (103)
This calculation is rather complex and is presented in detail in Appendix B.

The complete linearized operator is given in the following expression:

B ′ w i,j cijkl u k,l dv + B ′ w k,j σij u k,i dv + B ′ σij 1 3 (∇ x • u -∇ x • u)w i,j + 1 3 (∇ x • w -∇ x • w)u i,j + 1 9 (∇ x • u -∇ x • u)(∇ x • w -∇ x • w)δ ij + 1 9 (∇ x • u ∇ x • w -∇ x • u ∇ x • w)δ ij + 1 3 tr ∇ x u∇ x w -∇ x u∇ x w δ ij dv, (104)
where the notations introduced in Eqs. (B.5), (B.20) and (B.21) have been used. Note that, contrary to the F approach presented in de Souza Neto et al. [START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] De Souza Neto | F-bar based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF], the proposed method leads to a symmetric linearized operator.

Remark 1 From the implementational point of view, the same problems concerning the sparsity structure of the global system that were noted in the linear case arise here too. The "expected" contributions to the linearized operator, that is, the terms in Eq. ( 102), require a modification of the sparsity structure of the system equivalent to what is needed in the linear case. The additional "unexpected" terms require much more computational and storage effort because of terms involving projections of quantities that depend both on the weighting and solution spaces. However, the same strategy as in the linear case can be used. Although this requires the assembly of more "mixed" terms, the computational and storage cost of matrix-vector products is in our experience much smaller than the computation and assembly of the global matrix.

The discrete form of the consistent tangent operator given in Eq. ( 104) is presented in Appendix C.

Numerical examples in the nonlinear regime

The performance of the proposed F methodology is numerically assessed in this section. The first example, which is a nonlinear generalization of the Cook's membrane problem presented earlier, is widely used in the literature. Other two-and three-dimensional examples found in the literature are presented and results are compared with the other techniques mentioned in the introduction. 2 Cook's membrane. Evolution of the relative norm of the residual during the last Newton step for a 8 × 8 elements mesh with the proposed F technique with Q 1 /Q 0 .

Top right corner displacement / mm No F kref Q 3 No F kref Q 4 F Q 1 /Q 0 F kref Q 2 /Q 1 F kref Q 4 /Q 3 No F Q 4 No F Q 2 No F Q 1 No F Q 3 F Q 1 /Q 0 F Q 2 /Q 1 F kref Q 3 /Q 2 F Q 4 /Q 3 F Q 3 /Q 2

Number of elements per side

F kref Q 4 /Q 3 No F Q 1 No F kref Q 2 No F kref Q 3 No F kref Q 4 F Q 1 /Q 0 F kref Q 2 /Q 1 F kref Q 3 /Q 2 No F Q 2 F Q 1 /Q 0 No F Q 3 No F Q 4 F Q 4 /Q 3 F Q 3 /Q 2
We finally investigate the behavior of the consistently linearized operator described previously. We see in Table 2 the evolution of the relative Euclidian norm of the residual over the iterations for the last loading step for an 8 × 8 element mesh with the proposed F technique with Q 1 /Q 0 . This illustrates that the tangent matrix obtained from the consistent linearization of the problem results in quadratic convergence rate. The doubling of zeros behind the decimal point as we get very close to the solution is typically seen in our calculations.

Plane strain compression of a block

The next example is the plane strain compression of a rectangular block. This problem has been studied in Reference [START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF] as a severe benchmark for enhanced element formulations as they often exhibit locking in such situations [START_REF] Armero | On the formulation of enhanced strain finite elements in finite deformation[END_REF][START_REF] De Souza Neto | Remarks on the stability of enhanced strain elements in finite elasticity and elastocplasticity[END_REF]. The geometry, loading and boundary conditions are given in Figure 16. For symmetry reasons, only one half of the specimen needs to be considered. The material model is Neo-Hookean with a slightly different stored energy function than the one used in the previous example:

Ψ(J, C) = 1 2 µ (tr[C] -3) -µ ln J + 1 4 κ J 2 -1 -2 ln J (108)
Following Reference [START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF], the values for the material parameters are κ = 400889.806M P a and µ = 80.194M P a. The load factor p/p 0 is increased by steps of 0.5 with a reference value of p 0 = 20M P a (a unit thickness is considered). Note that a "dead" load is considered as shown in Figure 16, that is, one assumes a fixed vertical surface load specified on the reference configuration during the entire deformation of the structure. The Dirichlet boundary conditions are the following: no vertical displacement on the bottom and no horizontal displacement on the top. The quantity of interest is the compression level of the upper middle point. The evolution of the compression level versus the number of elements per edge for loading ratios from p/p 0 = 10 to p/p 0 = 60 is studied.

In Figure 17, we compare the results obtained with the proposed F approach with linear elements (i.e., Q 1 /Q 0 ) with the mixed element, referred to as Q 1 /P 0 (on the left), and with the enhanced element Q 1 SP , both from [START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF]. For all loading levels, convergence is attained for meshes of 16 × 16 elements or less. Note that, with coarse meshes, the F method performs better than both the Q 1 /P 0 and Q 1 SP elements at high compression levels while at low compression levels the situation is reversed. Results obtained with the proposed F approach for linear, quadratic and cubic NURBS (i.e., Q 1 /Q 0 , Q 2 /Q 1 and Q 3 /Q 2 ) are then compared in Figure 18. Although quadratic and cubic functions tend to produce a stiffer response than linears, we can see that in all cases we reach convergence for 16×16 element meshes. Quadratic and cubic NURBS produce better results for coarser meshes. Note that the jump from Q 2 /Q 1 to Q 3 /Q 2 only slightly improves the results as was observed for the previous example.

Three-dimensional compression of a block

The next example is the three-dimensional generalization of the previous one. It was studied in Reference [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF] as a benchmark for the three-dimensional version of the Q 1 SP element. This problem is of interest to test the performance of enhanced elements in three dimensional compression cases. The geometry, loading and boundary conditions are given in Figure 19, along with a typical plot of the Von Mises stress field on the final deformed configuration. For symmetry reasons, only one quarter of the specimen is considered.

The material model is Neo-Hookean with a somewhat different stored energy 40 17. Compression level in % versus number of elements per side for the two-dimensional plane strain block for various load levels. Left side: F Q 1 /Q 0 compared to the mixed u-p formulation Q 1 /P 0 (results from [START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF]). Right side: F Q 1 /Q 0 compared to the stabilized element Q 1 SP (results from [START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF]).

Q 1 SP F Q 1 /Q 0 Q 1 /P 0 p/p 0 =
function than the previous ones:

Ψ(J, C) = 1 2 µ (tr[C] -3) -µ ln J + 1 2 κ (ln J) 2 . ( 109
)
Following Reference [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF] the values for the material parameters are κ = 400889.806M P a and µ = 80.194M P a. Loading factors from p/p 0 = 20 to p/p 0 = 80 are studied with a reference value of p 0 = 4M P a. Note that a "dead" load is considered as shown in Figure 19, that is, one assumes a fixed vertical surface load specified on the reference configuration during the entire deformation of the structure. The Dirichlet boundary conditions are: no vertical displacement on the bottom surface and no horizontal displacement on the top surface. The quantity of interest is the compression level of the upper center point. The evolution of the compression level versus the number of elements per edge for various loading ratios is studied.

F Q 3 /Q 2 F Q 1 /Q 0 F Q 2 /Q 1 80 p/p 0 =
In Figure 20, we compare the results obtained with the proposed F approach with linear elements (i.e., Q 1 /Q 0 ) with the mixed element Q 1 /P 0 (on the left) and with the enhanced element Q 1 SP , both from [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF]. We can see that convergence is obtained for meshes of 8 × 8 × 8 elements and that with coarser meshes the results obtained with Q 1 /Q 0 are almost identical to the ones obtained with Q 1 /P 0 . Results obtained with linear, quadratic and cubic NURBS are presented in Figure 21. As for the two-dimensional case, quadratic and 20. Compression level in % versus number of elements per side for the three-dimensional block for various load levels. Left side: F

Q 1 SP F Q 1 /Q 0 Q 1 /P 0 compression level % Fig.
Q 1 /Q 0 compared to the mixed u-p formulation Q 1 /P 0 ([53]). Right side: F Q 1 /Q 0 compared to the stabilized ele- ment Q 1 SP ([53]).
cubic functions tend to give a stiffer response than linear functions. Nevertheless, results are sligthly better with the higher-order functions than with linear functions. The improvement obtained by going from quadratic to cubic NURBS is quite pronounced, especially for the coarse 2 × 2 × 2 element mesh.

Thick cylindrical shell

In the fourth example, we consider a geometry which cannot be represented exactly by piecewise linear polynomials. It involves the compression of a thick hollow cylinder, for which the exact initial geometry can be obtained with quadratic NURBS. This example was treated in [START_REF] Büchter | Three-dimensionnal extension of non-linear shell formulation based on the enhanced assumed strain concept[END_REF][START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF] and interest was focused on comparisons with shell formulations and shear locking. As we have seen in the Cooke's membrane problem, the proposed method is able to avoid both volumetric and shear locking. We note for this case that the usual Galerkin formulation using cubic and quartic functions also produced good results, without evidence of volumetric locking. The reason for this is that the traction-free cylindrical surfaces alleviate the tendency to lock. This tendency is also seen in plane stress analyses where there is no locking problem whatsoever. However, using the F method in conjunction with these high-order functions improves the results further. Therefore, we consider the thick cylinder with a compressible material in a bending case where shear locking is expected to occur to further assess this observation.

F Q 3 /Q 2 F Q 2 /Q 1 F Q 1 /Q 0
The geometry, boundary conditions, loading, and a typical result of the shear stress component σ yz of the Cauchy stress tensor on the final deformed configuration are given in Figure 22. As in [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF], we use only one quadratic element through the thickness, which is enough to capture the bending behavior. Because of symmetry, only one fourth of the structure is considered. The material model considered is the same as in the previous example with parameters κ = 240000M P a and µ = 6000M P a, which corresponds to a Poisson's ratio of ν = 0.4, therefore no volumetric locking is expected. Note that again a dead load is considered as shown in Figure 22, that is, one considers a fixed vertical line load specified on the reference configuration during the entire deformation of the structure. The quantity of interest is the vertical displacement of point A shown in Figure 22. 23. Vertical displacement of point A versus number of elements in the circumferential direction with and without F for various NURBS orders obtained from k-refinement, for Q 1 SP (results from [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF]) and for an enhanced assumed strain nonlinear shell element (results from [START_REF] Büchter | Three-dimensionnal extension of non-linear shell formulation based on the enhanced assumed strain concept[END_REF][START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF]).

F Q 2 /Q 1 No F Q 3 NL Shell No F Q 2 Q 1 SP F Q 3 /Q 2 18
F Q 2 /Q 1 No F kref Q 3 F kref Q 3 /Q 2 Q 1 SP NL Shell No F Q 2 Fig.
The results of the vertical displacement at point A are shown in Figure 23, where we compare the results obtained with and without F with C 1 -quadratic and C 2 -cubic NURBS, Q 1 SP and the enhanced assumed strain nonlinear shell element proposed in [START_REF] Büchter | Three-dimensionnal extension of non-linear shell formulation based on the enhanced assumed strain concept[END_REF], versus the number of elements in the circumferential direction. As in [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF], we used the following sequence of meshes (circumference × thickness × axis): 4×1×2, 8×1×4, 16×1×8, 32×1×16 that were obtained from k-refinement. We can see that all cases are converging to approximately the same value obtained in [START_REF] Reese | A new locking-free brick element technique for large deformation problems in elasticity[END_REF] with both Q 1 SP and shell elements, and that cubic functions give better results especially for the coarsest meshes. The use of F with quadratic functions improves the results by approximatively 50% for the coarse mesh, and 10% for the next mesh. Note also that, for cubic NURBS, the use of F improves the results for the coarse mesh by approximatively 10%, and that the two curves collapse for the third discretization.

Pinched torus

The last example again exploits the ability of NURBS to exactly represent conic sections. It consists of the pinching of a toroidal solid, and is similar to an example proposed in [START_REF] Chavan | Locking-free finite element methods for linear and nonlinear elasticity in 2D and 3D[END_REF]. The geometry, loading, boundary conditions and mesh are shown in Figure 24, and the material parameters are given in Table 3. The same Neo-Hookean material model as for the Cook's membrane is considered, which corresponds to a Poisson's ratio of ν = 0.4998. Due to symmetry condition, only one quarter of the structure is considered, with the corresponding symmetry boundary conditions applied. The quarter mesh with 4 × 16 × 2 elements (that is 4 elements in the "large" circumferential direction, 16 elements in the "small" circumferential direction and 2 elements in the radial direction) shown in Figure 24 with quadratic and cubic NURBS is used. The σ xx and σ zz components of the Cauchy stress tensor plotted on the final deformed configuration with and without F for C 1 -quadratic and C 2 -cubic functions are shown in Figure 25. The differences on the final deformed configurations show that Q 2 without F suffers form locking. Although F Q 2 /Q 1 and Q 3 without F look similar, we can see that F Q 3 /Q 2 improves the result. The stress contours for the component σ xx of the Cauchy stress tensor shows typical oscillations due to locking for Q 2 without F. These oscillations are not observed in the three other cases for this component. However, we can see on the bottom part of Figure 25, that oscillations are present for both quadratic and cubic meshes without F for the component σ zz of the Cauchy stress tensor. Note that, for both components and both orders of approximation, the results with F do not present such oscillations: the stresses are smooth and very similar for Q 2 /Q 1 and Q 3 /Q 2 .

The vertical displacement at point A for quadratic and cubic functions with and without F versus the loading step is shown in Figure 26. We also see that Q 2 suffers from locking and that at the final load the difference between Q 2 and Q 2 /Q 1 is approximatively of 50%. As we increase the order of approximation without F, we see that the results converge to those of Q 2 /Q 1 . Using C 2 cubic functions without F significantly improve the results which are very close to those of Q 2 /Q 1 . However, we can see that F Q 3 /Q 2 improves the result further. Note that, h-refining the mesh in the circumferential direction, that is, using a mesh of 8 × 16 × 2 elements, and F Q 2 /Q 1 allows us to obtain exactly the same results as for F Q 3 /Q 2 on the cruder mesh.

Elastic-plastic Cook's membrane

We present some preliminary results for elastic-plastic behavior at finite strain with the proposed F approach. We consider the same hyperelastic extension of J 2 -flow theory with poly-convex hyperelastic uncoupled stored energy function Fig. 28. Nonlinear elastic-plastic Cook's membrane with nonlinear isotropic hardening. Final equivalent plastic strain with F for a mesh of 16 × 16 elements using

σ xx /MPa F Q 2 /Q 1 no F Q 2 no F Q 3 F Q 3 /Q 2
Q 3 /Q 2 .
number of unknowns when increasing the order of approximation within the k -refinement process (contrary to p-refined finite elements) was also shown to be an interesting feature of isogeometric analysis.

This study demonstrates that the good properties of NURBS previously observed in other applications are also present with projection methods. The recent work on the use of high continuity functions (especially C 1 quadratics) compared with C 0 functions, in fluid mechanics [START_REF] Akkerman | The role of continuity in residual-based variational multiscale modeling of turbulence[END_REF][START_REF] Bazilevs | Isogeometric Analysis of turbulence and fluid-structure interaction[END_REF][START_REF] Bazilevs | Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[END_REF] and structural vibrations [START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF], gives us a better understanding of the effect of increased continuity on approximation. We anticipate that such properties could provide interesting improvements in large scale nonlinear solid mechanics.

As a first step, we proposed a new F projection method to treat the incompressibility constraint in the geometrically nonlinear regime. Although it shares features with the F method of de Souza Neto et al. [START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] De Souza Neto | F-bar based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF], we think of it as a more theoretically sound method since it is based on a modified minimum potential energy principle and leads to a symmetric consistent tangent which is not the case with the method of de Souza Neto et al. Moreover, the method does not make any assumption on the polynomial order of the approximation space and is therefore not limited to piecewise linear elements.

The numerical results show that the method compares well with previously published methods, especially for the Cook's membrane problem for which we obtained virtually identical results to those of [START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] De Souza Neto | F-bar based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF] with piecewise linear functions. We also investigated two-and three-dimensional compression prob-lems where the limits of enhanced strain elements are usually reached, and we observed good performance of the F isogeometric formulation even with low order functions. The ability of the method to overcome shear locking for quadratic and higher-order basis functions was also demonstrated although no particular attention was paid to that topic in the design of the method. Finally, some preliminary calculations considering an elastic-plastic behavior coupled with the proposed F method were presented, and show very promising results.

The proposed approach takes advantage of a simple pure displacement, projection method and high-order k -refined NURBS in the isogeometric analysis framework. It opens a new door to the use of precise CAD geometry representations and higher-continuity properties of NURBS-based isogeometric analysis in nonlinear structural mechanics applications. Areas of further research in structural mechanics needed include continuing the investigation of the approach initiated herein with elastic-plastic materials, and extending it to frictional contact. which can be introduced in Eq. (B.9), using the definition of the modified Cauchy stress tensor given in Eq. (96): SIJ U i,J W i,I = J σij w k,j u k,i .

(B.11)

The geometrical contribution to the linearized operator finally takes the expected form: Proceeding analogously, we can determine the next expression:

D(∇ x • W) • U = D π(J 1/3 ∇ x • W) J 1/3 • U. (B.19)
Following Eq. (B.5), we introduce the following notations:

∇ x • w ∇ x • u = π J 1/3 ∇ x • w ∇ x • u J 1/
D(∇ x • W) • U = 1 3 ∇ x • w ∇ x • u -tr [∇ x w∇ x u] - 1 3 ∇ x • u ∇ x • w. (B.22)
The corresponding discrete internal forces can be expressed as: (J 1/3 N A ,L F -1 Li , ÑC ) M -1 CD ( ÑD , J -1/3 SIJ FkI F kJ )w A i .

Fig. 1 .

 1 Fig. 1. Basis functions of order 0, 1, 2 for an open, uniform knot vector Ξ = {0, 1, 2, 3, 4, ...}.

( a )Fig. 3 .

 a3 Fig. 3. Piecewise quadratic B-Spline curve in R d . a) Control point locations are denoted by 's. b) Knots, which define the mesh by partitioning the curve into elements, are denoted by 's.

Fig. 4 .

 4 Fig.[START_REF] Bazilevs | Isogeometric Analysis of turbulence and fluid-structure interaction[END_REF]. Curve, control points, knot vectors and basis functions after knot insertion and order elevation of a quadratic curve. Control points are denoted by 's.

Fig. 6 .

 6 Fig.6. Basis functions Q p /Q p-1 , p = 1, 2, 3, for a one-dimensional patch of four elements. All cases attain continuity C p-1 /C p-2 .

Fig. 7 .

 7 Fig.7. Basis functions Q p /Q p-1 , p = 2, 3, with a point of C 0 continuity within the patch, for a one-dimensional patch of four elements. All cases attain continuity C p-1 /C p-2 , except at the repeated knot ξ = 0.75 where the continuity is only C 0 .

Fig. 8 .

 8 Fig.8. Geometry, loading, boundary conditions, material parameters and quantity of interest for the plane strain Cook's membrane.

Fig. 10 .

 10 Fig. 10. Cook's membrane. Vertical displacement of top right corner versus number of degrees of freedom with and without B for various NURBS orders obtained from k-refinement. The case "No B Q1" is offscale.

Fig. 11 .

 11 Fig. 11. Geometry, loading, boundary conditions and material parameters for the infinite plate with a hole under tension.

Fig. 12 .

 12 Fig. 12. Infinite plate with a hole: sequence of meshes produced by h-refinement.

Q 2 Fig. 14 .F Q 2

 2142 Fig. 14. Vertical displacement of top right corner versus number of elements per edge with and without F for various NURBS orders obtained from k-refinement.

Fig. 15 .

 15 Fig. 15. Vertical displacement of top right corner versus number of degrees of freedom with and without F for various NURBS orders obtained from k-refinement. Note that Q 1 results without F are not shown on the figure.

Fig. 16 .

 16 Fig. 16. 2D plane strain compression test: geometry, boundary conditions, loading, and typical results of Von Mises Stress field on the deformed configuration.

20 p

 20 Fig.[START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF]. Compression level in % versus number of elements per side for the two-dimensional plane strain block for various load levels. Left side: F Q 1 /Q 0 compared to the mixed u-p formulation Q 1 /P 0 (results from[START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF]). Right side: F Q 1 /Q 0 compared to the stabilized element Q 1 SP (results from[START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF]).

Fig. 18 .

 18 Fig. 18. Compression level in % versus number of elements per side for the two-dimensional plane strain block for various load levels. Comparison of F results for various NURBS orders obtained from k-refinement.

Fig. 19 .

 19 Fig. 19. 3D compression test: geometry, boundary conditions, loading, and typical results of the Von Mises stress field on the deformed configuration.

Fig. 21 .Fig. 22 .

 2122 Fig. 21. Compression level in % versus number of elements per side for the three-dimensional block for various load levels. Comparison of F results for various NURBS orders obtained from k-refinement.

  in circumferential direction Vertical displacement at point A / mm

Fig. 24 .

 24 Fig. 24. Pinched torus: geometry, quarter mesh, loading and boundary conditions.

  torus: material properties and boundary conditions.

Fig. 25 .

 25 Fig. 25. Pinched torus: component (a) σ xx and (b) σ zz of the Cauchy stress tensor on the deformed configuration with and without F for C 1 -quadratic functions and C 2 -cubic functions (note that the scales have been reduced to see the interesting features).

  j σij u k,i dv. (B.12)To determine the second term in Eq. (B.9), we need to compute:DW i,I • U, (B.13)whereW i,I = 1 3 (∇ x • W -∇ x • W) FiI FiJ + α FiJ W i,I (B.14)Rewriting the divergence of W as∇ x • W = W i,i = W i,I F -1 Ii = tr(∇ X W F -1 ), (B.15) we get D(∇ x • W) • U = tr(∇ X W • D(F -1 )) U. (B.[START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF]) Since we have the following equation for the directional derivative of the inverse of the deformation gradient:D(F -1 ) • U = -F -1 ∇ X U F -1 , (B.17)Eq. (B.16) can be simplified as:D(∇ x • W) • U = -tr [∇ x w∇ x u] . (B.[START_REF] Cottrell | Isogeometric analysis of structural vibrations[END_REF]

  .., ξ n+p+1 }, let ξ ∈ [ξ k , ξ k+1 [ be a desired new knot. The new n + 1 basis functions are formed recursively, using Eqs. (2) and (3), with the new knot vector Ξ = {ξ 1 , ..., ξ k , ξ, ξ k+1 , ..., ξ n+p+1 }. The new n + 1 control points {B 1 , ..., B n+1 } are formed from the original control points {B 1 , ..., B n }, by

Table 1 Comparison

 1

	Finite Element Analysis	Isogeometric Analysis
	Nodal points	Control points
	Nodal variables	Control variables
	Mesh	Knots
	Basis interpolates nodal	Basis does not interpolate
	points and variables	control points and variables
	Approximate geometry	"Exact" geometry
	Polynomial basis	NURBS basis
	Gibbs phenomena	Variation diminishing
	Subdomains	Patches
	Compact support
	Partition of Unity
	Isoparametric concept
	Affine covariance
	Satisfy patch tests

of Finite Element and Isogeometric Analysis properties.

also called stored energy or strain energy function

Acknowledgements T. Elguedj was supported by the French Délégation Générale pour l'Armement and the Centre National pour la Recherche Scientifique. Y. Bazilevs was partially supported by the J.T. Oden ICES Postdoctoral Fellowship at the Institute for Computational Engineering and Sciences. Y. Bazilevs, V.M. Calo and T.J.R. Hughes were partially supported by the Office of Naval Research under Contract No. N00014-03-0263, Dr. Luise Couchman, contract monitor, and Sandia National Laboratories under Contract No. 114166. This support is gratefully acknowledged.

Mesh parameter h Relative error in L 2 norm of stress

Mesh parameter h Mesh parameter h Relative error in energy norm [START_REF] Büchter | Three-dimensionnal extension of non-linear shell formulation based on the enhanced assumed strain concept[END_REF]. Plate with a hole. Convergence curves in relative L 2 norm of displacement, relative L 2 norm of stress and relative energy norm with and without B for various NURBS order obtained from k-refinement.

Cook's membrane

The setup of the nonlinear Cooke's membrane example is similar to the one previously presented in the geometrically linear case. The geometry, loading and boundary conditions are given in Figure 8. This example has been solved in References [START_REF] Armero | On the formulation of enhanced strain finite elements in finite deformation[END_REF][START_REF] Chavan | Locking-free finite element methods for linear and nonlinear elasticity in 2D and 3D[END_REF][START_REF] Cesar De Sa | New enhanced-strain elements for incompressible problems[END_REF][START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF][START_REF] De Souza Neto | F-bar based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking[END_REF][START_REF] Kasper | A mixed-enhanced strain method. Part II: geometrically nonlinear problems[END_REF][START_REF] Klass | A stabilized mixed finite element method for finite elasticity. Formulation for linear displacement and pressure interpolation[END_REF][START_REF] Maniatty | Higher order stabilized finite element method for hyperelastic finite deformation[END_REF][START_REF] Reese | A new stabilization technique for finite elements in non-linear elasticity[END_REF][START_REF] Simo | Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes[END_REF] with enhanced assumed strain methods, the F method of de Souza Neto et al., and stabilized mixed formulations.

The material model used is Neo-Hookean following the additive decomposition of the stored energy function given in Eq. ( 77). The isochoric and volumetric parts of Ψ are (see, e.g., [START_REF] Simo | Computational Inelasticity[END_REF]):

where κ is the bulk modulus and µ the shear modulus. Following Reference [START_REF] De Souza Neto | Design of simple low order finite elements for large strain analysis of nearly incompressible solids[END_REF], the values for the material parameters are κ = 40.0942 × 10 4 M P a and µ = 80.1938M P a which corresponds to a Poisson ratio of ν = 0.4998. The load value is identical to the one used in the linear case, that is, F = 100N/mm. The quantity of interest is again the vertical displacement of the top right corner of the plate. The convergence as a function of the number of elements per edge and total number of degrees of freedom are used as criteria to investigate the performance of the method. The results for the vertical displacement versus the number of elements per edge are given in Figure 14 and the ones for the vertical displacement versus the total number of degrees of freedom are given in Figure 15.

The results are similar to the ones obtained with the B method in the linear case. The converged result obtained for quartic NURBS functions is very close to values obtained with previously published results. Without any treatment of the incompressibility condition, increasing the polynomial order of the NURBS functions improves the results, however, we can see that C 1 quadratic NURBS still suffer from locking for coarse meshes: we need a 32 × 32 elements mesh to attain the converged result. The use of F significantly improves the results, and good convergence is attained even with Q 1 /Q 0 . Note that the results obtained with Q 1 /Q 0 are better than Q 2 and that Q 2 /Q 1 give better results than Q 3 .

For high-order NURBS with F, we obtain converged results with very coarse meshes. The results in term of displacement versus total number of degrees of freedom show again the advantages of k -refinement. Note that in terms of unknowns, no and nonlinear isotropic hardening as in Simo [START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part I: continuum formulation[END_REF][START_REF] Simo | A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part II: computational aspects[END_REF] and Simo and Hughes [START_REF] Simo | Computational Inelasticity[END_REF]. It consists of a Neo-Hookean model for the elastic part, and an associative flow rule based on a Von Mises yield criterion with isotropic hardening following a saturation law for the plastic part. The nonlinear isotropic hardening rule is defined by the following equation:

where σ 0 is the initial flow stress, σ ∞ is the saturation flow stress, δ is the saturation exponent, K is the linear hardening coefficient, and α is the equivalent plastic strain.

The Neo-Hookean material model and boundary conditions are identical to the ones previously used for this problem, and the load value is chosen to be F = 5N/mm. The material parameters are those used in Simo and Armero [START_REF] Simo | Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes[END_REF], and are presented in Table 4.

The results, shown in Figure 27 for the vertical displacement of the top right point versus number of elements per side, compare well to those obtained in [START_REF] Armero | On the formulation of enhanced strain finite elements in finite deformation[END_REF][START_REF] Simo | Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes[END_REF]. Similar convergence behavior, when increasing the order of approximation with the proposed F method, to those obtained in the nearly incompressible linear and nonlinear elastic cases is obtained. One can also see the equivalent plastic strain for a 16 × 16 element mesh with F using cubic NURBS functions in Figure 28.

Conclusions

In this paper we have presented B and F projection methods for nearly incompressible small deformation and large deformation elasticity in the framework of isogeometric analysis. These methods make use of the higher continuity property of high-order NURBS obtained from k -refinement, a unique feature of isogeometric analysis. We first presented an extension of the well known B method in this context. This allowed us to concentrate on the choice of projection spaces and techniques to solve problems with higher-order functions. The results obtained showed the efficacity of the projection methods to solve nearly incompressible problems with a pure displacement formulation using higher-order functions with higher continuity. The small increase in the

A Directional derivative of the modified Green-Lagrange strain tensor

To obtain the directional derivative of Ē, we first need to compute the derivative of the modified deformation gradient F:

Using the linearity of the projection operator and expanding the former expression, we get:

For simplicity, we make use of the notation introduced in Eq. ( 75):

and therefore Eq. (A.2) becomes:

The definition of the modified Green-Lagrange strain tensor in terms of the modified deformation gradient is used in conjunction with the previous equation:

which finally allows us to write:

B Consistent linearized operator

We begin by recalling Eq. ( 103):

We focus on the first term of this equation. We define the modified fourth order material tensor as we previously defined Eq. ( 91), the modified second Piola-Kirchhoff stress tensor:

Using the result from Appendix A (i.e., Eq. (A.6)), the first term of Eq. (B.1) becomes CIJKL FiJ W i,I FjL U j,K . (B.3) As done previously, we can push forward this equation to the current configuration using:

where we have introduced for brevity the following notation in both the reference and current configuration:

The push forward of Eq. (B.3) is:

We can consider the classical push forward on the current configuration of the modified material tensor in terms of "bar" quantities, that is:

and finally obtain the expected material contribution to the linearized operator on the current configuration:

We now consider the second term in Eq. (B.1), which can in turn be split into two terms:

The first term represents the geometrical contribution to the linearized operator in terms of "bar" quantities, and the second term is the one responsible for the additional contribution. We first focus on the geometrical contribution.

Using the result obtained in Appendix A (i.e., Eq. (A.6)), we have:

We can gather all terms to obtain the following expression for Eq. (B.13):

Using this expression we finally obtain the last unexpected contribution to the linearized operator:

C Discrete consistent linearized operator

To express the discrete version of the consistent linearized operator given in Eq. (104), we start by introducing the following:

The "material" contribution to the discrete consistent linearized operator (i.e., the first term in Eq. ( 104)) is

The "geometrical" contribution to the discrete consistent linearized operator (i.e., the second term in Eq. ( 104)), and parts of the "unexpected" contributions are:

Finally, the remaining terms of the "unexpected" contribution are:

Li , ÑC ) M -1 CD ( ÑD , J -1/3 SIJ FkI F kJ )u B j + 1 9 w A i (J 1/3 N A ,K F -1 Ki N B ,L F -1 Lj , ÑC ) M -1 CD ( ÑD , J -1/3 SIJ FkI F kJ )u B j .