exact geometry

Abstract The concept of isogeometric analysis is proposed. Basis functions generated from NURBS (Non-Uniform Rational B-Splines) are employed to construct an exact geometric model. For purposes of analysis, the basis is refined and/or its order elevated without changing the geometry or its parameterization. Analogues of finite element h - and p -refinement schemes are presented and a new, more efficient, higher-order concept, k -refinement, is introduced. Refinements are easily implemented and exact geometry is maintained at all levels without the necessity of subsequent communication with a CAD (Computer Aided Design) description. In the context of structural mechanics, it is established that the basis functions are complete with respect to affine transformations, meaning that all rigid body motions and constant strain states are exactly represented. Standard patch tests are likewise satisfied. Numerical examples exhibit optimal rates of convergence for linear elasticity problems and convergence to thin elastic shell solutions. A k -refinement strategy is shown to converge toward monotone solutions for advection–diffusion processes with sharp internal and boundary layers, a very surprising result. It is argued that isogeometric analysis is a viable alternative to standard, polynomial-based, finite element analysis and possesses several advantages.


Introduction
In this paper we introduce a new method for the analysis of problems governed by partial differential equations such as, for example, solids, structures and fluids. The method has many features in common with the finite element method and some features in common with meshless methods. However, it is more geometrically based and takes inspiration from Computer Aided Design (CAD). A primary goal is to be geometrically exact no matter how coarse the discretization. Another goal is to simplify mesh refinement by eliminating the need for communication with the CAD geometry once the initial mesh is constructed. Yet another goal is to more tightly weave the mesh generation process within CAD. In this work we introduce ideas in pursuit of these goals.
It is interesting to note that finite element analysis in engineering had its origins in the 1950s and 1960s. Aerospace engineering was the focal point of activity during this time. By the late 1960s the first commercial computer programs (ASKA, NASTRAN, Stardyne, etc.) appeared. Subsequently, the finite element method spread to other engineering and scientific disciplines, and now its use is widespread and many commercial programs are available. Despite the fact that geometry is the underpinning of analysis, CAD, as we know it today, had its origins later, in the 1970s and 1980s. A highly-recommended introductory book, with historical insights, is Rogers [1]. This perhaps explains why the geometric representations in finite element analysis and CAD are so different. Major finite element programs were technically mature long before modern CAD was widely adopted. Presently, CAD is a much bigger industry than analysis. Analysis is usually referred to as Computer Aided Engineering (CAE) in market research. It is difficult to precisely quantify the size of the CAE and CAD industries but current estimates are that CAE is in the $1-$2 billion range and CAD is in the $5-$10 billion range. The typical situation in engineering practice is that designs are encapsulated in CAD systems and meshes are generated from CAD data. This amounts to adopting a totally different geometric description for analysis and one that is only approximate. In some instances mesh generation can be done automatically but in most circumstances it can be done at best semi-automatically. There are still situations in major industries in which drawings are made of CAD designs and meshes are built from them. It is estimated that about 80% of overall analysis time is devoted to mesh generation in the automotive, aerospace, and ship building industries. In the automotive industry, a mesh for an entire vehicle takes about four months to create. Design changes are made on a daily basis, limiting the utility of analysis in design if new meshes cannot be generated within that time frame. Once a mesh is constructed, refinement requires communication with the CAD system during each refinement iteration. This link is often unavailable, which perhaps explains why adaptive refinement is still primarily an academic endeavor rather than an industrial technology.
The geometric approximation inherent in the mesh can lead to accuracy problems. One example of this is in thin shell analysis, which is notoriously sensitive to geometric imperfections; see Fig. 1. The sensitivity to imperfections is shown in Fig. 1b in which the buckling load of a geometrically perfect cylindrical shell is compared with shells in which geometric imperfections are introduced with magnitudes of 1%, 10%, and 50% of the thickness. As may be seen, there is a very considerable reduction in buckling load with increased imperfection.
Sensitivity to geometry has also been noted in fluid mechanics. Spurious entropy layers about aerodynamic shapes were the bane of compressible Euler solvers in the 1980s and 1990s. The problem and its solution were identified in the thesis of Barth [2]. Piecewise linear approximations of geometry were the root cause. Smooth geometry completely eliminated the entropy layers even when the flow fields were approximated by linear elements on the curved geometry; see Fig. 2. This result explains why methods which employ smooth geometric mappings are widely used in airfoil analysis (see [3]). It is also well known in computational fluid dynamics that good quality boundary layer meshes significantly improve the accuracy of computed wall quantities, such as pressure, friction coefficient, and heat flux; see Fig. 3.
The construction of finite element geometry (i.e., the mesh) is costly, time consuming and creates inaccuracies. It is clear from the smaller size of the CAE industry compared with the CAD industry that the most fruitful direction would be to attempt to change, or replace, finite element analysis with something more CAD-like. This direction was taken in the development of the RASNA program Mechanica, in which exact geometry in conjunction with a p-adaptive finite element procedure was utilized. However, the lack of satisfaction of the isoparametric concept led to theoretical questions which were addressed in later versions of the code by abandoning the exact geometry in favor of high-order polynomial approximations [7]. The use of a fixed polynomial approximation to geometry has been shown by Szabó et al. [8] to be limiting. As solution polynomial order is increased, the error plateaus at some level and cannot be further reduced (see Fig. 4). The seriousness of this result is compounded by the fact that computed quantities defined on boundaries are usually the most important ones in engineering applications, and this is where geometric errors are most harmful. Furthermore, most finite element analyses are still performed with low-order elements for which geometric errors are largest. The success of RASNA, which was later acquired by Parametric Technology Corporation (PTC), a CAD company, was due to its tight linkage with CAD  [4] and (b) buckling of cylindrical shell with random geometric imperfections [5]. solution space is piecewise linear, while geometry is piecewise quadratic. Smooth geometry avoids spurious entropy layers associated with piecewise-linear geometric approximations (from [6]).
geometry and, perhaps more importantly, its consequent ability to provide adaptive p-refinement and thus more reliable results. The present methodology is similarly inspired, but attempts to more faithfully adhere to CAD geometry and eliminate the finite element polynomial description entirely. (The p-method is described in Szabó and Babuška [9] and Szabó et al. [8].) The approach we have developed is based on NURBS (Non-Uniform Rational B-Splines), a standard technology employed in CAD systems. We propose to match the exact CAD geometry by NURBS surfaces, then construct a coarse mesh of ''NURBS elements''. These would be solid elements in three-dimensions that exactly represent the geometry. 1 This is obviously not a trivial task and one that deserves much study but when it can be accomplished it opens a door to powerful applications. Subsequent refinement does not require any further communication with the CAD system and is so simple that it may facilitate more widespread adoption of this technology in industry. There are analogues of h-, p-, and hp-refinement strategies, and a new, higher-order methodology emerges, k-refinement, which seems to have advantages of efficiency and robustness over traditional p-refinement. All subsequent meshes retain exact geometry. Throughout, the isoparametric philosophy is invoked, that is, the solution space for dependent variables is represented in terms of the same functions which represent the geometry. For this reason, we have dubbed the methodology isogeometric analysis. NURBS are not a requisite ingredient in isogeometric analysis. We might envision developing isogeometric procedures based on ''A-patches'' (see [11][12][13][14][15][16]) or ''subdivision surfaces'' (see [17][18][19]). However, NURBS seem to be the most thoroughly developed CAD technology and the one in most widespread use.
The body of this paper begins with a tutorial on B-splines (B-splines are the progenitors of NURBS), followed by one on NURBS. We then describe an analysis framework based on NURBS. This is followed by sample applications in linear solid and structural mechanics and some introductory calculations in fluids, namely, ones involving classical test cases for the advection-diffusion equation. Various refinement strategies are studied and, in cases for which exact elasticity solutions are available, optimal rates of convergence are attained. The structural problems include some applications to thin shells modeled as solids. The approach is seen to handle these situations remarkably well. In the fluid calculations, we employ the SUPG formulation and consider difficult test cases involving internal and boundary layers. We observe that, by employing highorder, k-refinement strategies, convergence toward monotone solutions is obtained. This surprising result seems to contradict numerical analysis intuitions and suggests the possibility of linear difference methods that are simultaneously robust and highly accurate. We close with conclusions and suggestions for future work.

Knot vectors
NURBS are built from B-splines. The B-spline parametric space is local to ''patches'' rather than elements. Patches play the role of subdomains within which element types and material models are assumed to be uniform. A knot vector in one dimension is a set of coordinates in the parametric space, written N ={n 1 , n 2 , ..., n n+p+1 }, where n i 2 R is the ith it knot, i is the knot index, i = 1,2, ..., n + p +1,p is the polynomial order, and n is the number of basis functions which comprise the B-spline.
Remark. The convention we will adopt is that the order p = 0, 1, 2, 3, etc., refers to constant, linear, quadratic, cubic, etc., piecewise polynomials, respectively. This is the usual terminology in the finite element literature. What we refer to as ''order'' is usually referred to as ''degree'' in the computational geometry literature.
If knots are equally-spaced in the parametric space, they are said to be uniform. If they are unequally spaced, they are non-uniform. More than one knot can be located at the same coordinate in the parametric space. These are referred to as repeated knots. A knot vector is said to be open if its first and last knots appear p + 1 times. Open knot vectors are standard in the CAD literature. In one dimension, basis functions formed from open knot vectors are interpolatory at the ends of the parametric space interval, [n 1 , n n+p+1 ], and at the corners of patches in multiple dimensions but they are not, in general, interpolatory at interior knots. This is a distinguishing feature between knots and ''nodes'' in finite element analysis.

Basis functions
B-spline basis functions are defined recursively starting with piecewise constants (p =0) For p = 1,2,3, ..., they are defined by Derivatives with respect to spatial coordinates may be computed by way of standard techniques described in Hughes [20,Chapter 3]. An initial example of the results of applying (1) and (2) to a uniform knot vector is presented in Fig. 5. Note that, for p = 0 and 1, the basis functions are the same as for standard piecewise constant and linear finite element functions, respectively. However, for p P 2, they are different. Quadratic B-spline basis functions (and NURBS basis functions, as will be shown later) are identical but shifted. This is in contrast with quadratic finite element functions which are different for internal and end nodes. This ''homogeneous'' pattern continues as we go to higher-order B-splines and may result in significant advantages in equation solving over finite element functions, which are quite ''heterogeneous''.
An example of quadratic basis functions for an open, non-uniform knot vector is presented in Fig. 6. Note that the basis functions are interpolatory at the ends of the interval and also at n = 4, the location of a repeated knot, where only C 0 -continuity is attained. Elsewhere, the functions are C 1 -continuous. In general, basis functions of order p have p À 1 continuous derivatives. If a knot is repeated k times, then the number of continuous derivatives decreases by k. When the multiplicity of a knot is exactly p, the basis function is interpolatory. (1) They constitute a partition of unity, that is, "n X n i¼1 N i;p ðnÞ¼1: ð3Þ (2) The support of each N i,p is compact and contained in the interval [n i , n i+p+1 ].
(3) Each basis function is non-negative, that is, N i,p (n) P 0,"n. Consequently, all of the coefficients of a mass matrix computed from a B-spline basis are greater than, or equal to, zero.

B-spline curves
B-spline curves in R d are constructed by taking a linear combination of B-spline basis functions. The coefficients of the basis functions are referred to as control points. These are somewhat analogous to nodal coordinates in finite element analysis. Piecewise linear interpolation of the control points gives the so-called control polygon. In general, control points are not interpolated by B-spline curves. Given n basis functions, N i,p , i =1,2,..., n, and corresponding control points B i 2 R d ; i ¼ 1; 2; ...; n, a piecewise-polynomial Bspline curve is given by An example is shown in Fig. 7 for the quadratic basis functions considered previously. Note that the curve is interpolatory at the first and last control points, due to the fact that the knot vector is open, and also at the sixth control point, due to the fact that the multiplicity of the knot n = 4 is equal to the polynomial order. Note also that the curve is tangent to the control polygon at the first, last, and sixth control points. The curve is C pÀ1 = C 1 -continuous everywhere except at the location of the repeated knot, n = 4, where it is C pÀ2 = C 0 -continuous. (1) They have continuous derivatives of order p À 1 in the absence of repeated knots or control points.
(2) Repeating a knot or control point k times decreases the number of continuous derivatives by k.
(3) An affine transformation of a B-spline curve is obtained by applying the transformation to the control points. 2 We refer to this property as affine covariance.

h-refinement: knot insertion
The analogue of h-refinement is knot insertion. Knots may be inserted without changing a curve geometrically or parametrically. Given a knot vector N ={n 1 ,n 2 , ..., n n+p+1 }, let n 2½n k ; n kþ1 ½ be a desired new knot. The new n + 1 basis functions are formed recursively, using (1) and (2), with the new knot vector N ¼fn 1 ; n 2 ; ...; n k ; n; n kþ1 ; ...; n nþpþ1 g. The new n + 1 control points, fB 1 ; B 2 ; ...; B nþ1 g, are formed from the original control points, {B 1 , B 2 , ..., B n }, by where Knot values already present in the knot vector may be repeated as above but as described in Section 2.2, the continuity of the basis will be reduced. Continuity of the curve is preserved by choosing the control points as in (5) and (6). Each unique internal knot value may appear no more than p times or the curve becomes discontinuous.
An example of knot insertion is presented in Fig. 8. The original curve consists of quadratic B-splines. The knot vector is N = {0, 0, 0, 1, 1, 1}. The curve is shown on the left with basis functions below. A new knot is inserted at n ¼ 0:5. The new curve, shown on the right, is geometrically and parametrically identical to the original curve, but the basis functions, below the curve, and control points are changed. There is one more of each. This process may be repeated to enrich the solution space by adding more basis functions of the same order while leaving the curve unchanged. This subdivision strategy is seen to be analogous to the classical h-refinement strategy in finite element analysis.

p-refinement: order elevation
The polynomial order of basis functions may be increased without changing the geometry or parameterization. It is important to note that each unique knot value in N must be repeated to preserve discontinuities in the pth derivative of the curve being elevated. The number of new control points depends on the multiplicities of existing knots. This strategy of order elevation is an analogue of p-refinement in finite element analysis.
As is the case of h-refinement by way of knot insertion, the solution space spanned by the order elevated basis functions contains the space spanned by the original functions. Thus, it is possible to order elevate without changing the geometry of the B-spline curve. Less obviously, it can be done so as to leave the parameterization of the curve in tact. The process for doing this involves subdividing the curve into many Bézier curves by knot insertion (see [1] or [21] for a discussion of Bézier curves), order elevating each of these individual segments, and then removing the unnecessary knots to combine the segments into one, order-elevated, B-spline curve. Several efficient algorithms exist which combine the steps so as to minimize the computational cost of the process. We omit the details for the sake of brevity. For a thorough treatment, see Piegl and Tiller [21]. (Note that the CAD community refers to this process as ''degree elevation''.) An example of order elevation is depicted in Fig. 9. The original curve and quadratic basis functions, shown on the left, are the same as considered in the previous example. This time the multiplicity of the knots is increased by one. The numbers of control points and basis functions each increase by one. The locations of the control points change, but the elevated curve is geometrically and parametrically identical to the original curve. There are now four cubic basis functions. The locations of control points for this elevated curve are different than those in the previous example (cf. Fig. 8).

k-refinement
An alternative order elevation strategy takes advantage of the fact that the processes of knot insertion and order elevation do not commute. If a unique knot value, n, is inserted between two distinct knots in a curve of order p, the number of continuous derivatives of the basis functions at n is p À 1. As described above, if we subsequently elevate the order to q, the multiplicity of every distinct knot value (including the knot just inserted) is increased so that discontinuities in the pth derivative of the basis are preserved, that is, the basis still has p À 1 continuous derivatives at n. If instead we elevated the order of the original curve to q and only then inserted a unique knot value, the basis would have q À 1 continuous derivatives at n. We refer to this latter procedure as k-refinement. It has no analogue in standard finite element analysis. We believe the concept of k-refinement is important and potentially a superior approach to high-precision analysis than p-refinement. 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. In addition, there is a proliferation in the number of nodes because C 0 -continuity is maintained in the refinement process. In k-refinement, there is a homogeneous structure within patches and growth in the number of control variables is limited. This will be made clear with an example but first we need a definition. Let us define an ''element'' in one dimension as the span between two distinct knot values. The number of elements in a curve will then be the number of non-zero knot spans in the knot vector (e.g., the domain associated with the knot vector N = {0, 0, 0, 1, 2, 3, 3, 4, 4, 4} consists of four elements). This definition anticipates our NURBS-based finite element implementation, described in Section 3.
With this notion of an element, consider the refinement processes depicted in Fig. 10. Assume the initial domain consists of one element and p + 1 basis functions (assuming an open knot vector), which we then hrefine until we have n À p elements and n basis functions. We then perform order elevation by p-refinement, maintaining continuity at the p À 1 level. This requires replicating each distinct knot value, adding a basis function in each element and so increasing the total number of basis functions by n À p to 2n À p. After a total of r order elevations of this type, we have (r +1)n À rp basis functions, where p is still the order of our original basis functions. This is seen to be a large number of functions when one considers that in most cases of practical interest the number of elements will be quite a bit larger than the order of the basis. By comparison, begin with the same one element domain and proceed by k-refinement. That is, order elevate r times adding only one basis function at each refinement, then h-refine until we have n À p elements as before. The final number of basis functions is n + r, each having r + p À 1 continuity. This amounts to an enormous savings as n + r is considerably smaller than (r +1)n À rp. Keep in mind too that in d dimensions these numbers are raised to the d power. There may also be other advantages in that smoother derivatives may lead to more accurate physical quantities, such as strains and stresses, and increased smoothness may, surprisingly, lead to better capturing of thin layers. An example of this phenomenon is presented later on. Of course, if the physical situation dictates a certain lower level of continuity at a knot value ,1 ,1,1},p=2 = {0, 0, 0, ,1,1,1},p=2 . k-refinement takes advantage of the fact that knot insertion and order elevation do not commute. (a) Base case of one linear element. (b) Classic p-refinement approach: knot insertion followed by order elevation results in seven piecewise quadratic basis functions that are C 0 at internal knots. (c) New k-refinement approach: order elevation followed by knot insertion results in five piecewise quadratic basis functions that are C 1 at internal knots.
(e.g., the corners in the geometry, discontinuous material properties, etc.) this can always be incorporated into the process by knot duplication, so no generality is lost.

B-spline surfaces
Given a control net {B i,j }, i = 1,2, ..., n, j =1,2,..., m, and knot vectors N ={n 1 , n 2 , ..., n n+p+1 }, and H ¼fg 1 ; g 2 ; ...; g mþqþ1 g, a tensor product B-spline surface is defined by where N i,p and M j,q are basis functions of B-spline curves. For purposes of numerically integrating arrays constructed from B-splines, ''elements'' are taken to be knot spans, namely, [n i , n i+1 ] · [g j , g j+1 ]. See Fig. 11 for an illustration of a standard bi-unit parent element and its image in physical space. Integrals are pulled back to the parent element by the classical change-of-variables formula and standard Gaussian quadrature rules are employed; see [20, Chapter 3].

Rational B-splines
Desired geometric entities in R d can be obtained by projective transformations 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 quadratic curves. This is illustrated in Fig. 12 in which a circle in R 2 is constructed from a piecewise quadratic B-spline curve in R 3 . The projective transformation of a B-spline curve yields a rational polynomial of the form C R (n)=f(n)/g(n), where f and g are piecewise polynomials. The construction of a rational B-spline curve in R d proceeds as follows. Let fB w i g be a set of control points for a B-spline curve in R dþ1 with knot vector N. These are referred to as the ''projective'' control points for the desired NURBS curve in R d . The control points in R d are derived from the projective control points by the following relations: where (B i ) j is the jth component of the vector B i , etc. and w i is referred to as the ith weight.InFig. 12a, the weights are the vertical coordinates of the control points defining the piecewise quadratic B-spline curve in R 3 . The rational basis functions and NURBS curve are given by Rational surfaces and solids are defined analogously in terms of the rational basis functions As an example, Fig. 13 shows a control net and the corresponding NURBS surface description of a torus. The paraphernalia describing the constructions in Figs. 12 and 13 are provided in Appendices A.1 and A.2, respectively. Important properties of NURBS are: (1) NURBS basis functions form a partition of unity.
(2) The continuity and support of NURBS basis functions are the same as for B-splines.
(3) Affine transformations in physical space are obtained by applying the transformation to the control points, that is, NURBS possess the property of affine covariance.

NURBS as a basis for analysis
An analysis framework based on NURBS consists of the following items and features: (1) A mesh for a NURBS patch is defined by the product of knot vectors. For example, in three dimensions, a mesh is given by N Â H Â Z. (4) The control points associated with the basis functions define the geometry. (5) The isoparametric concept is invoked, that is, the fields in question (e.g., displacement, velocity, temperature, etc.) are represented in terms of the same basis functions as the geometry. The coefficients of the basis functions are the degrees-of-freedom, or control variables. (6) Mesh refinement strategies are developed from a combination of knot insertion and order elevation techniques. These enable analogues of classical h-refinement and p-refinement methods, and the new possibility of k-refinement. (7) Arrays constructed from isoparametric NURBS patches can be assembled into global arrays in the same way as finite elements; see Hughes [20,Chapter 2]. Compatibility of NURBS patches is attained by employing the same NURBS edge and surface representations on both sides of patch interfaces; see Fig. 14a. This gives rise to a standard continuous Galerkin method. Refinement necessarily propagates from patch to patch. There are two alternatives corresponding to Fig. 14b. One is the discontinuous Galerkin method, which can be employed at the patch level. Compatibility is enforced weakly by the variational formulation. The other is to utilize constraint equations for the control points and control variables to attain pointwise compatibility at patch interfaces. See Kagan et al. [22] for a procedure of this kind in the context of B-splines. Local refinement is an important and challenging research topic. (8) The easiest way to set Dirichlet boundary conditions is to apply them to the control variables. In the case of homogeneous Dirichlet conditions, this results in exact, pointwise satisfaction. In the case of inhomogeneous Dirichlet conditions, the boundary values must be approximated by functions lying within the NURBS space. This amounts to ''strong'' but approximate satisfaction of the boundary con- Remark. It is well known that typical finite element interpolation functions oscillate in attempting to fit discontinuous data. An example is illustrated in Fig. 15a where Lagrange polynomials of orders three, five, and seven interpolate a discontinuity represented by eight data points in R 2 . Note that as the order is increased, the amplitude of the oscillations increases. This is sometimes referred to as Gibbs phenomena. NURBS behave very differently when the data are viewed as control points. In Fig. 15b the NURBS curves are monotone, illustrating the variation diminishing property of NURBS (see [1,Chapter 3]). This property has advantages in representing sharp layers. An example illustrating this will be presented later.  A summary of similar and dissimilar finite element and isogeometric analysis concepts is presented in Table 1. A salient feature of isogeometric analysis, shared by meshless methods, is the non-interpolatory nature of the basis. Other shared features are the partition of unity property and compact support. For a state-of-the-art summary of meshless methods, see [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].

Isogeometric structural analysis
Isoparametric NURBS patches represent all rigid body motions and constant strain states exactly. This follows from the properties of NURBS with respect to affine transformations. Structures assembled from compatible NURBS patches pass standard ''patch tests''. (Patch tests are described in Chapters 3 and 4 of Hughes [20]. They are considered by engineers to be practical ways of assessing whether or not finite elements are convergent and whether or not they have been programmed correctly.) We have verified this assertion in numerous test cases (not shown). In what follows, we present numerical solutions for linear elastic solids and structures. The Galerkin formulation of linear elasticity is employed (see [20,Chapter 2]). The first example is two-dimensional and the remainder are three-dimensional. Several of the calculations involve thin shells. These are modeled as NURBS solids and no shell assumptions are employed. Gaussian quadrature is used on elements. The rule of thumb that we typically employed is to use the lowest-order rule that would exactly integrate the integrand assuming the NURBS are B-splines of the same polynomial order and the Jacobian determinants are constants. This represents an approximation in the case of NURBS. To assess its validity, we performed tests in which we systematically increase the number of quadrature points. For sufficiently fine meshes no differences in results were discernible. However, coarse meshes required more integration points due to large variations in the geometrical mapping. More research needs to be done to determine a robust strategy covering all situations.
Remark. We employed both direct and iterative linear algebraic equation solvers. The direct solver had a profile architecture and was typically more efficient for shell structures than the iterative, diagonally preconditioned, conjugate gradient solver. However, this was not the case for NURBS solids. The code we developed was restricted to a single processor and many of the finer-mesh cases could not be solved with the direct solver because they exceeded available memory. The aspect ratios of elements in some cases were quite large and the number of equations approached a quarter of a million. Nevertheless, in all cases the iterative procedure converged without difficulty. We did not test the corresponding finite element cases but we would have been surprised if iteration behaved this reliably in those cases too. Consequently, we conjecture that the NURBS cases may be better conditioned than the corresponding finite element cases. The evidence we have for this is rather skimpy but we think it warrants further investigation.

Infinite plate with circular hole under constant in-plane tension in the x-direction
In this two-dimensional example, we present in some detail the NURBS analysis of a problem in solid mechanics having an exact solution. We also systematically explore h-and k-refinement. The infinite plate is modeled by a finite quarter plate. The exact solution [41, pp. 120-123], evaluated at the boundary of the finite quarter plate, is applied as a Neumann boundary condition and is given here for reference: where T x is the magnitude of the applied stress for the infinite plate case. The setup is illustrated in Fig. 16.
R is the radius of the hole, L is the length of the finite quarter plate, E is YoungÕs modulus, and m is Pois-sonÕs ratio. A rational quadratic basis is the minimum order capable of representing a circular hole. The coarsest mesh, N Â H, is defined by the knot vectors Exact traction 16. Elastic plate with a circular hole: problem definition.  The exact geometry is represented with only two elements, as shown in Fig. 17a. The corresponding control net is shown in Fig. 17b. A summary of the geometric data is given in Appendix A.3. A repeated control point is responsible for the upper-left hand corner. Fig. 18 depicts three different views of the mesh and support of two of the basis functions. These are the index space, parameter space, and physical space views. The index space illustrates how only parts of the basis functions are supported by the physical space, identified by the dark-bordered rectangle (i.e., [n 3 , n 5 ] · [g 3 , g 4 ]) in Fig. 18a. This is due to the open knot vectors. Consequently, knot spans near the boundary have zero measure and may be ignored in the element formation phase.
Contour plots of the same two basis functions are shown in Fig. 19. The first six meshes used in the analysis are shown in Fig. 20. The scheme used to select knots for insertion is described in Appendix B.
Contour plots of results obtained on meshes 1, 4, and 7 are presented in Fig. 21. The applied stress is T x = 10 and the contours show that the stress concentration of r xx = 30 at the edge of the hole (i.e., at r = R, h = 3/2p) is obtained as the mesh is refined.
Convergence results in the L 2 -norm of stresses are shown in Fig. 22. The cubic and quartic NURBS are obtained by order elevation of the quadratic NURBS on the coarsest mesh. Since the parameterization of  the geometrical mapping does not change, the h-refinement algorithm generates identical meshes for all polynomial orders. As a result, the continuity of the basis is C pÀ1 everywhere, except along the line which joins the center of the circular edge with the upper left-hand corner of the domain (see Fig. 17a). There it is C 1 as is dictated by the coarsest mesh employing rational quadratic parameterization. The mesh parameter, h max , is defined as the maximum distance, in physical space, between diagonally opposite knot locations. As  can be seen, the L 2 -convergence rates of stress for quadratic, cubic, and quartic NURBS are approximately 2, 3, and 4, respectively. This is the best one could reasonably hope for. The convergence of the standard deviation of weights for rational quadratic elements is depicted in Fig.  23. As may be seen, the standard deviation converges to zero, meaning the weights converge to a constant within each element and the NURBS approach polynomial B-splines. The geometric way to see this is to note that the mesh in physical space (i.e., R 2 ) can be ''lifted'' to projective space (i.e., R 3 ). As the elements in projective space shrink in size, the projective control points move closer together. As the weights are just the third component of the projective control points, within an element they necessarily converge to a constant value. This geometric argument suggests the rate of convergence should be O(h max ), as observed in Fig. 23.

Solid circular cylinder subjected to internal pressure loading
The problem specification is shown in Fig. 24. Plane strain conditions are assumed to hold in the axial direction. It is a simple matter to obtain an exact solution to this problem, assuming the pressure varies at most circumferentially (see [41,[117][118][119]). We provide the exact solution for the case with constant pressure as a reference: Meshes developed from h-refinement are shown in Fig. 25. The NURBS paraphernalia for the coarsest mesh is given in Appendix A.4.
Results for the constant pressure case, employing quadratic NURBS, are presented in Fig. 26.I nFig. 26a the axisymmetric deformation of the cylinder is clearly apparent for Mesh 4. Axisymmetric response is attained for all meshes. Errors in radial displacement are plotted in Fig. 26b. For the first (coarsest) mesh, the maximum error through the wall thickness is approximately 1%. For the second mesh it is approximately 0.1% and for the third 0.01%. Order elevated solutions employing cubic and quartic NURBS are, for all practical purposes, exact on all meshes. Hence, results for these cases are not presented. In order to determine convergence rates, a somewhat more complex loading was considered in which the internal pressure was assumed to vary as cos 2h. Quadratic, cubic, and quartic NURBS were considered. The cubic and quartic cases were obtained from the quadratic case by k-refinement, in which case the degree of continuity was increased to C 2 and C 3 , respectively. The rates of convergence of the error measured in the energy norm are presented in Fig. 27. As in the example of the plate with a circular hole, the rates of convergence for quadratic, cubic, and quartic NURBS elements are 2, 3, and 4, respectively.

Solid ''horseshoe'' subjected to equal and opposite in-plane flat edge displacements
The problem setup is presented in Fig. 28a. The top surfaces are displaced in the directions shown. This introduces an asymmetry to the loading that is visible in the stress contours. Meshes 1, 3, and 5 are depicted  in Figs. 28b and 29. The standard h-refinement procedure, as described in Appendix B, is employed to obtain finer meshes. It is clear for this case that good quality, refined meshes can be simply generated, once an exact parameterized geometry is constructed in the form of the first mesh. Fig. 30 shows sample plots of stress contours for a quadratic NURBS solution on Mesh 5. Note the smoothness of the stress contours.

Thin cylindrical shell with fixed ends subjected to constant internal pressure
The problem setup and a radial displacement profile are shown in Fig. 31. Note that the fixed ends create boundary layers which are difficult to accurately capture with low order finite element methods. The exact shell theory solution is given in Timoshenko and Woinowsky-Krieger [42, pp. 476-477], for plane stress and is provided as a reference below  uðxÞ¼À PR 2 Et ð1 À C 1 sin bx sinh bx À C 2 cos bx cosh bxÞ x 2ðÀL=2; L=2Þ; ð23Þ C 1 ¼ sin a cosh a À cos a sinh a sinh a cosh a þ sin a cos a ; ð24Þ  The YoungÕs modulus and PoissonÕs ratio in this solution need to be replaced by E 1Àm 2 and m 1Àm , respectively, to account for the fixed-end conditions assumed here. The geometry of the shell is shown in Fig. 32. Note that the radius to thickness ratio is 100. The meshes are depicted in Figs. 33 and 34. The first four surface meshes are shown in Fig. 33. Note the added refinement in the region of the boundary layer. We wish to emphasize that, despite the shell being very thin, we are modeling it with solid elements. (For an excellent comprehensive review of approaches to shell modeling, see Bischoff et al. [43].) The through-thickness mesh resolution for Mesh 1 is shown in Fig. 34. Note that there are two elements in the radial direction and four in the circumferential direction. As the surface mesh is refined, the number of elements in the circumferential direction increases accordingly. However, the number of elements in the radial direction is fixed at two throughout the refinement process. The functions employed for this problem are quadratic NURBS. Their appearance in the radial, or through-thickness, direction is presented in Fig.  34a. Radial displacement contours are presented in Fig. 35 for Meshes 2-5. Note that, for all meshes, perfectly axisymmetric response is obtained. Note also the appearance of boundary layers. The convergence of the radial displacement profile is shown on Fig. 36. Mesh 1 is too coarse to represent both the boundary layers and the plateau between. Mesh 3 picks up the plateau but the boundary layers are still not accurately     captured. The Mesh 5 solution is indistinguishable from the exact shell theory solution. In the detail on the right, the exact shell solution and Mesh 5 solution are seen to overlap in the boundary layer region.

Shell obstacle course
The so-called shell obstacle course consists of three problems, the Scordelis-Lo roof, the pinched hemisphere, and the pinched cylinder. These problems, and their relevance to the assessment of shell analysis procedures, have been discussed extensively in the literature. The problem descriptions presented in Figs. 37-39, are adapted from Felippa [44] and Belytschko et al. [45]. As in the previous example, two quadratic  NURBS elements are employed in the through-thickness direction (see Fig. 34), whereas h-refinement and k-refinement are utilized for surface meshing. Quadratic through quintic surface NURBS are employed in the convergence analysis of all cases. In one case, the pinched hemisphere, a one-element surface solution, starting with rational quadratics, the lowest-order NURBS capable of exactly representing spherical geometry, and culminating with 10th-order NURBS, is used to assess convergence. This analysis has the flavor of what are usually referred to as ''spectral methods,'' which are accurate and efficient procedures, typically utilized for performing detailed studies of geometrically simple but physically complex phenomena, such as turbulence (see [46,47] for detailed descriptions and applications of spectral methods). Convergence is assessed by comparing the displacement of certain points in the shell with benchmark solutions presented in Belytschko et al. [45]. Sample contour plots of the solutions for quadratic elements and the finest meshes studied are presented in Figs. 40-42. Note that in each case the contours are very smooth.

Scordelis-Lo roof
The Scordelis-Lo roof is subjected to gravity loading. The ends are supported by fixed diaphragms and the side edges are free (see Fig. 37). The vertical displacement of the midpoint of the side edge is the quantity used to assess convergence. The second, fourth, and sixth meshes used in the study are shown in Fig. 43. These meshes have 2, 8, and 32 surface elements per side, respectively. Due to symmetry, only one quadrant  is meshed. Convergence of the displacement to the benchmark value is shown in Fig. 44. In all cases, convergence is quite rapid. For the higher-order cases, namely, quartic and quintic, even one element provides a very accurate solution.

Pinched hemisphere
In the pinched hemisphere, equal and opposite concentrated forces are applied at antipodal points of the equator. The equator is otherwise considered to be free (see Fig. 38). The second, fourth, and sixth meshes are shown in Fig. 45. Due to symmetry, only one quadrant is meshed. Convergence of the displacement under the inward directed load is presented in Fig. 46. The quadratic case converges very slowly, which is not surprising as quadratic, fully-integrated, solid finite elements are known to ''lock'' in shell analysis. Cubic solid finite elements also exhibit locking in similar circumstances but in the present case cubic NURBS behave reasonably well. Fig. 47 presents convergence of the displacement for one surface element meshes. Notice that the lowest-order meshes lock but eventually accurate results are obtained. One 10thorder NURBS surface element is seen to provide an essentially exact result. To assess whether there is any tendency to oscillate, displacement in the direction of the inward directed point load is plotted for the single 10th-order NURBS surface element case in Fig. 48. As is evident, the displacements are very smooth and monotone.

Pinched cylinder
The pinched cylinder is subjected to equal and opposite concentrated forces at its midspan (see Fig. 39). The ends are supported by rigid diaphragms. This constraint results in highly localized deformation under the loads (see Fig. 42). Only one octant of the cylinder is used in the calculation due to symmetry. The second, fourth, and sixth meshes are shown in Fig. 49. Convergence of the displacement under the load is    presented in Fig. 50. The NURBS elements converge to a very slightly softer solution than the benchmark solution. This may be due to transverse shear effects. It is well known that, as long as the characteristic surface element dimension is large compared with the thickness, formulations which permit transverse shear deformations typically closely approximate formulations which satisfy the Kirchhoff constraint (i.e., zero transverse shear). When this trend reverses, that is as the surface element dimension approaches zero, holding the thickness constant, the displacement under a concentrated load grows, and converges to infinity (see [48] for elaboration).
Remark. We believe that two quadratic NURBS through the thickness are unnecessary for typical thin shell analysis. We performed some tests with one quadratic NURBS through the thickness and the results were indistinguishable when compared with the two-quadratic-NURBS case. However, when we reduced to linear variation through the thickness, convergence to correct solutions was not obtained. Classic shell theory hypotheses, such as invoking the plane stress condition in the through-thickness direction, may be sufficient to correct the deficiency of linear through-thickness displacement variation. Such a formulation would certainly be competitive with traditional shell element formulations which employ displacement and rotation degrees of freedom at a reference surface. Using only displacement degrees of freedom (i.e., control variables) in the NURBS case considerably simplifies shell analysis, especially in non-linear analysis wherein rotations are no longer vectorial and additive, but require a multiplicative group structure. A further simplification might be just to use a NURBS surface and employ ''rotationless'' formulations, such as the one for plates described in Engel et al. [49]. In this reference a discontinuous Galerkin formulation is proposed but all interface discontinuity (i.e., ''jump'') terms in it disappear if C 1 , or higher, continuity is satisfied. This is simply attained with NURBS but very difficult to achieve in finite element analysis. Rotationless shell elements have recently gained popularity in computational mechanics (see [17][18][19][49][50][51][52][53]). The formulation of Engel et al. [49] has also been proposed as being an appropriate basis for so-called ''strain gradient theories''. These theories also require C 1 -continuity, and NURBS would appear to be naturally suited to them.

Hemispherical shell with a stiffener
We consider the hemispherical shell with stiffener presented in Rank et al. [10] who solved the problem using a finite element method and p-refinement strategy. The geometry and data are depicted in Fig. 51. The structure is subjected to gravity loading and external pressure, and vertical displacements are set to zero on the bottom surface of the stiffener. Only a quarter of the domain is modeled due to symmetry. The initial mesh is constructed using rational quadratic NURBS and is shown in Fig. 52. Isotropic p-refinement was performed to assess convergence. Polynomial orders 2-6 were employed. Figs. 53 and 54 show the vertical displacement and von Mises stress on the deformed configuration. The von Misses stress is given by Fig. 51. Hemispherical shell with stiffener. Problem description from Rank et al. [10].
where r vm is the von Mises stress, r ij is the stress tensor and r 0 ij is the deviatoric stress. The results in Figs. 53 and 54 came from a quadratic NURBS simulation. There is qualitative agreement with the results given in Rank et al. [10]. The Euclidean norm of the displacement and the von Mises stress were calculated at points A-D, identified in Fig. 51. Since points A-D are located on the boundary of the computational domain, values of r ij imposed by the boundary conditions were used in (27), where applicable, as opposed to the computed quantities. Results for p =2 to p = 6 are plotted in Figs. 55-62 along with results from Rank et al. [10] for their finest simulation, that is, p = 8, which we take as a reference. Good agreement in the converged displacements is observed. Figs. 56 and 57 indicate that at points B and C, the displacements for the NURBS discretization converge from the ''stiff'' side (consistent with the results of the previous section). Agreement in converged values of the von Mises stress is also apparent.

Isogeometric fluids analysis
We wish to perform preliminary assessments of isogeometric analysis in fluids. The fluid equations are fundamentally different than the structural equations in that advection dominated flow phenomena are characterized by essentially skew-symmetric differential operators, the antithesis of the symmetric, definite operators which characterize typical structural analysis models. Another complicated flow feature is sharp layers involving very strong gradients. Within these layers, diffusive behavior prevails and the equations transition to being predominantly symmetric-definite. As a result, successful fluid analysis formulations need to automatically account for the local competition between advective and diffusive effects, and adjust the discretization accordingly. The standard Galerkin formulation optimally treats symmetric-definite operators but produces unstable discretizations of skew-symmetric operators. This is well known in finite element analysis and, unfortunately, carries over to the NURBS-based approach. The first finite element methodology exhibiting uniform stability and convergence behavior across the full range of advective and diffusive phenomena was SUPG, a ''stabilized method'' which had its origins in the late 1970s and early 1980s. (A good summary of the very early literature is contained in Brooks and Hughes [54]. Recent, stateof-the-art literature on stabilized methods is presented in [55][56][57][58][59][60][61][62][63][64][65][66].) SUPG is a residual-based modification of the Galerkin method which increases stability without degrading accuracy. It is sufficient for most flow phenomena but it is not a monotone technology because its higher-order accuracy normally precludes monotone behavior except in sufficiently resolved cases. There is a well known ''barrier theorem'' that asserts that the class of higher-order accurate, linear difference methods cannot be monotone, and SUPG falls into that class. Consequently, SUPG by itself is not a sufficiently robust technology for shock wave phenomena in fluids but in combination with so-called ''discontinuity-capturing'' operators has proved to be an industrial-strength technology for applications with shocks (see [67]). SUPG possesses a rich mathematical theory and there is an extensive literature on its mathematical properties beginning with Johnson et al. [68]; see Hughes et al. [69] and references therein for a review of the literature. The identification of stabilized methods as multiscale methods was made in Hughes [70] and Hughes et al. [71]. These works also explored the link between stabilized methods and subgrid-scale modeling [69], an important direction in current research.
In this section we investigate the ability of the isogeometric approach, in conjunction with SUPG, to solve challenging test cases for the advection-diffusion equation. Isogeometric analysis is fundamentally a higher-order approach and one would not expect good behavior in situations with unresolved interior and boundary layers. In fact, the Gibbs phenomena noted for polynomial-based finite element methods tends to become more pronounced as polynomial order is increased. This is the reason that most practical fluids formulations employ lower-order, typically constant and linear, interpolation of flow variables. However, the variation diminishing property of Dirichlet boundary condition specification, plus the notion of krefinement, leads to some remarkable results in the case of NURBS. This is illustrated by the first example.

Advection skew to the mesh with outflow Dirichlet boundary conditions
The problem setup is described in Fig. 63. Here, a is the advective velocity magnitude, j is the diffusivity, and L is the side length of the domain. The Peclet number, Pe = aL/j =10 6 . When this number is greater than one, advection dominates and diffusion is only important in very small layers. In the present case, diffusion is important in a region of thickness O(Pe À1 ln Pe) in the outflow boundary layers and O(Pe À1/2 ln Pe) in the internal layer (see [72, p. 468]). In all calculations the mesh is uniform, consisting of a 20 · 20 grid of square elements, with element side length h = 1/20 = 0.05. Refinement is performed by the k-method, and solutions from p =1 to p = 12 are calculated. In all cases, the standard SUPG formulation is used with s = h a /2a, where h a is the element length in the direction of the flow velocity which, in the present case, is simply, h a = h/max{cos h, sin h}.    The boundary condition is set by specifying the control variables. On the top and right edges of the domain, all control variables are set to 0 and the boundary condition is exactly satisfied along these edges. On the bottom, the control variable corresponding to the lower right-hand corner is set to 0 and the remainder are set to 1. The result is that the boundary value is identically 1 up to the last element in which it smoothly decreases to 0 at the corner. The left-hand-side boundary is more interesting. If we think of our control variables as control points in R 3 Fig. 65. For p = 1, the boundary condition is interpolated, whereas for p > 1 it is fit to the control variables in monotone fashion as the variation diminishing property of B-splines prevents the curve from over-and under-shooting. We wish to emphasize that k-refinement produces nonnested solution spaces, which prevents us from having the exact same boundary condition at each stage of the refinement process. As a result of this technique, the discontinuity ''smears'' about the location 0.2. For odd polynomial orders with p 6 9, a control point falls directly on 0.2, whereas for even polynomial orders one does not (see Fig. 64). For this reason we have plotted the even and odd cases separately. When p is larger than 9, the ''fringing'' due to the open knot vectors disrupts this pattern. Clearly, a better scheme for setting the boundary control variables is required. One possibility is to obtain the desired boundary condition by setting up a projection problem over the boundary, for which a suitable norm needs to be chosen. Our next computational example will demonstrate that technique. Nevertheless, the main point of the present study is to assess the ability of NURBS (in this case, B-splines because of the simplicity of the domain) to deal with unresolved boundary and interior layers. We first consider the case in which h =45°. The results for p =1top = 12 are presented in Figs. 66-68. Two views are presented for each p, one in which the plotting routine sampled the solution with a 100 · 100 grid of uniformly distributed points and one in which it is sampled with a 21 · 21 uniform grid. In the former case the plot is Phong shaded and in the latter, it is represented by bilinear interpolation on each element and the element edges are drawn. The philosophy behind the dual views is that the 100 · 100 grid plots are a more faithful rendering of the higher-order cases, whereas the 21 · 21 point piecewise bilinear interpolates are the type of plots that have appeared in numerous research articles over the years and these may be more easily visually compared with results in the literature. The p =1,h =45°, case is unusual in that the selection of s, the SUPG parameter, is the tuned value that captures the boundary layer very well. In general, the boundary layer is very poorly captured with p =1    unless s takes on a specific value. As one examines the results, it is clear that they improve as p increases and are converging toward monotone results with quite sharp layers. One would perhaps expect that oscillations would increase with increasing p but this is not the case. This is certainly due in part to the smearing of the boundary condition but we suspect that the high continuity of the basis obtained through k-refinement plays a part in this as well. Further studies need to be undertaken to clarify these issues.
It is interesting to view the control variables (see Fig. 69). The solution for p = 8 is nearly monotone, whereas the corresponding control variables are quite oscillatory. It is not clear what physical information is provided by the control variables. However, as the mesh is h-refined, the control net will approach the solution surface, that is, the basis becomes interpolatory in the limit as h ! 0. The 20 · 20 mesh is obviously far from this limit.
The h = tan À1 (2) case (i.e., h % 63.4°) is presented in Figs. 70-72. This problem represents possibly the most severe standard test for the advection-diffusion equation. Refinement is performed by the k-method, and solutions from p =1top = 12 are calculated. This time, for p = 1, the overshoot in the outflow boundary layer is approximately 45% but, as the order is elevated by k-refinement, convergence toward monotone results occurs. This angle of flow seems to be somewhat more difficult than the h =45°case. Nevertheless, the same conclusion may be drawn: Sufficiently high-order k-refinement produces extraordinary results in dealing with sharp internal and boundary layers. The maximum percentage overshoot vs. polynomial order for both cases is plotted in Fig. 73. The convergence toward monotone solutions is apparent in both cases.

Advection of a sine hill in a rotating flow field
The problem statement is shown in Fig. 74. The flow is circular about the center of a square domain [À1, 1] · [À1, 1] with velocity components given by a 1 = Àx 2 and a 2 = x 1 . Diffusivity is taken to be 10 À6 , which makes the problem advection-dominated. Homogeneous essential boundary conditions are imposed around the perimeter of the domain, while the u = Àsin (2px) condition is imposed on the ''slit,'' as shown in Fig. 74. Since the B-spline space does not contain trigonometric functions, H 1 0 -minimization was used to set the control point values on the ''slit''. Computations were done on a 30 · 30 uniform mesh with p = 2 and p = 6. Continuity was maintained at the C 1 and C 5 levels, respectively, everywhere in the domain except along the lines x = 0 and y = 0 where C 0 continuity was maintained. This choice was made to facilitate imposition of boundary conditions. Results are presented in Fig. 75. The exact solution is essentially a pure advection of the boundary condition along the circular streamlines. Both p = 2 and p = 6 produce excellent results.

Conclusions
We have introduced a new analysis framework, called isogeometric analysis, which is based on NURBS. It has many similarities with finite element analysis, but also important differences. A key feature is to represent geometry exactly by NURBS elements and then invoke the isoparametric concept to define field variables, such as displacement, temperature, etc. The coarsest mesh encapsulates the exact geometry. This means that mesh refinement is simply accomplished by reindexing the parametric space. The refinement process can proceed without interaction with the CAD system, a distinct advantage over finite element procedures. The fact that finite element mesh refinement strategies require interaction with the CAD system at each stage may be the reason that adaptive mesh refinement, despite its benefits, is still primarily an academic activity and one that has not significantly penetrated the industrial sector. There are NURBS analogues of finite element h-refinement and p-refinement, and there is also a variant of p-refinement, which we call k-refinement, in which the continuity of functions is systematically increased. This seems to have no analogue in traditional finite element analysis but is a feature shared by some meshless methods. The approach is fundamentally higher-order. For example, in order to represent circles, cylinders and spheres, rational polynomials of at least quadratic order are necessary. The generation of refined NURBS bases of all orders is facilitated by simple recursion relationships. The versatility and power of recursive NURBS basis representations are truly remarkable. Equation systems generated by NURBS tend to be more homogeneous than those generated by higher-order finite elements and this may have some benefit in equation solving strategies. NURBS satisfy a ''variation diminishing'' property. For example, they give monotone fits to discontinuous control data and become smoother as order is increased, unlike Lagrange interpolatory polynomials which oscillate more violently as order is increased. NURBS of all orders are non-negative pointwise. This means that every entry of the NURBS mass matrix is non-negative. These properties are not attained in finite element analysis. On the other hand, NURBS are not interpolatory. They are fit to nets of control points and control variables. This aspect is less transparent to deal with than the corresponding finite element concepts of interpolated nodal points and nodal variables but somewhat similar to the situation for meshless methods.
We performed initiatory calculations with isogeometric analysis procedures on some linear structural and fluids problems. h-, p-, and k-refinement strategies were investigated. In the structural applications, convergence to exact solutions was shown to occur at optimal rates for two-and three-dimensional solids. We performed shell analysis with three-dimensional solid NURBS elements. Convergence to thin shell solutions and benchmark results occurred in all cases. The ability to accurately resolve shell boundary layers appears to be particularly noteworthy.
We applied a stabilized SUPG formulation of isogeometric analysis to standard test cases involving the advection-diffusion equation in order to assess its potential for fluid dynamics. A stringent two-dimensional test case, advection skew to a mesh with outflow Dirichlet boundary conditions, was studied on a fixed mesh using k-refinement. As order was increased, solutions improved and converged toward monotone results with sharp layers. This surprising result seems to contradict numerical analysis preconceptions about high-order accurate methods.
There seem to be numerous opportunities for research in isogeometric analysis. Some areas of particular interest are: (1) integrating isogeometric analysis with CAD systems; (2) developing isogeometric mesh generators based on state-of-the-art technologies such as CUBIT (see [73]); (3) developing triangular and tetrahedral isogeometric elements and associated unstructured meshing techniques (see [74]); (4) developing procedures which account for trimmed surface description (this is a significant challenge, the approach of Höllig [75] may be useful, see Fig. 76); (5) developing isogeometric surface elements for traditional and rotationless shell analysis [17][18][19][49][50][51][52][53] and boundary integral formulations, and combining isogeometric surface descriptions with meshless methods (exact surface descriptions would seem to be particularly easy to generate from CAD descriptions); (6) applying isogeometric technology to thin plate and shell theories, strain-gradient theories, and other higher-order theories; (7) developing hierarchical isogeometric bases for multiscale analysis; (8) applying isogeometric analysis to fluid mechanics and, in particular, boundary layer turbulence (the pioneering studies [76][77][78][79] should be mentioned in this regard); (9) applying isogeometric elements to problems of engineering interest; (10) developing a mathematical theory of convergence and error analysis (we have already made some headway on this topic and out initial results will be reported upon in a forthcoming publication); (11) applying isogeometric analysis to contact problems and, in particular, frictional contact and sliding (see Fig. 77); (12) developing isogeometric dynamic analysis, including time-harmonic and explicit, implicit, and space-time time-stepping approaches (see Fig. 78); (13) applying isogeometric analysis to shape optimization (see Fig. 79); (14) generalizing to other important application areas, such as, for example, acoustics and electromagnetics; (15) various non-linear applications, etc.
In conclusion, we believe the isogeometric approach has considerable potential in practical problem solving and is a promising alternative to current analysis procedures.

Acknowledgement
This research was supported by Office of Naval Research Contract N00014-03-0263, Dr. Luise Couchman, contract monitor. J.A. Cottrell was supported by a VIGRE Fellowship. This support is gratefully acknowledged. We also wish to thank Ernst Rank and Alexander Dü ster for access to their unpublished data, papers, and helpful comments, and Marc Halpern for insight into the CAD and CAE industries.

A.2. Toroidal surface
The toroidal surface shown in Fig. 13b is created by sweeping a circular NURBS curve of radius 1, as defined in the previous section, around a circular path of radius 5, also defined by a NURBS curve. In this instance, the g coordinate describes our initial circular curve and the n coordinate traverses the path along which we sweep. The basis is rational quadratic in each of the two directions, with knot vectors given by N ¼f0; 0; 0; 1; 1; 2; 2; 3; 3; 4; 4; 4gð A:2Þ and H ¼f0; 0; 0; 1; 1; 2; 2; 3; 3; 4; 4; 4g: ðA:3Þ As N is the same as in our previous example, the set of functions, {N i,2 (n)}, are the same as in Table A.1. As N ¼ H, the set {M j,2 (g)} have the same form, but with independent variable g rather than n.
The control net for the torus may be thought of as given by ''sweeping'' the control polygon for a circle of radius 1 about the control polygon for a circle of radius 5. See Fig. 13a. The resulting set of control points is given in Table A.3, with weights in Table A.4. The rational surface basis functions are formed from the one-dimensional basis functions and weights using (13). A.3. Plate with circular hole The two-element coarse mesh shown in Fig. 17a and defined by the control net in Fig. 17b has an interesting feature. The circle requires using a rational quadratic basis, but this in turn necessitates introducing a discontinuity in the first derivative to create the opposite corner. This could be done either by repeating knots and thus creating a basis that is only C 0 along the line between the initial two elements, or by placing two control points at the same location in physical space. We chose the latter option in order to ensure that our basis had C 1 continuity throughout the interior of the domain. The only place where the basis is not C 1 is on the boundary itself, at the location of the repeated control points. As stated previously, this mesh uses knot vectors N ¼f0; 0; 0; 0:5; 1; 1; 1gð A:4Þ N 1,2 (n)( 2 n À 1) 2 0 N 2,2 (n)2 n(2À3n)2 ( n À 1) 2 N 3,2 (n)2 n 2 2(3n À 1)(1 À n) N 4,2 (n)0 ( 2 n À 1) 2  and H ¼f0; 0; 0; 1; 1; 1g: ðA:5Þ The corresponding one-dimensional basis functions are given in Tables A.5 and A.6. The control net and weights are given in Tables A.7 and A.8, respectively.

A.4. Cylindrical solid
The cylinder shown in Fig. 25 has an inner radius of r in = 1, and an outer radius of r out = 2. The length of the cylinder is L = 5. The coarse mesh uses quadratic basis functions in all three parametric directions. The n coordinate traverses the circumferential direction, the g coordinate traverses the thickness, and the f coordinate traverses the length. As each of these knot vectors have appeared in a previous example, analytic expressions for the corresponding one-dimensional basis functions may be found in the above tables. The control net and weights are given in Table A.9. Rational solid basis functions are defined by combining the weights and one-dimensional basis functions using (14).
A hollow cylindrical solid can be developed with a single NURBS patch by mapping opposite faces of the parametric domain to the same physical location in the obvious way. For open knot vectors, the control points describing the two faces are coincident. This leads to a C 0 -continuous solution across the interface.

Appendix B. Knot selection during h-refinement
The difficulty encountered in h-refinement is how to determine where to insert knots. Knot locations may always be chosen by the user to increase the resolution in the vicinity of specific features but an automated global refinement strategy is desired. The parameterization of a NURBS patch in physical space is dependent not only on the locations of control points and the polynomial order of the basis functions but also on the weights. In practice, placing a new knot in the middle of an existing knot span, thus splitting the element into two equal parts in the parametric domain, often does not even come close to splitting the element in half in physical space. Complicating matters further is the fact that knot insertion is inherently a global process. For example, consider a mesh in R 2 defined by N Â H. Inserting a knot n 2½n i ; n iþ1 into N may split the image of the element edge [n i , n i+1 ] · g 1 in physical space in half, but fail to split the image of [n i , n i+1 ] · g m+q+1 in half, and so no unique knot value will exactly halve every element in the row. Still, in meshes described by one patch, refinement requires selecting one value. The approach we have developed has worked well in the examples we have considered, and is described below for the two-dimensional case.
When inserting a knot into [n i , n i+1 ], perform the following steps: (1) Find the physical coordinates corresponding to (n i , g 1 ) and (n i+1 , g 1 ), and calculate the midpoint of the chord connecting them; call this point P 1 . Likewise, find P 2 , the midpoint of the chord from (n i , g m+q+1 ) and (n i+1 , g m+q+1 ). See Fig. B.1.
(2) Perform point inversion to find n 1 and n 2 , the n coordinates minimizing the distance between P 1 and P 2 and points on the images of [n i , n i+1 ] · g 1 and [n i , n i+1 ] · g m+q+1 , respectively. See Fig. B.1. This is done using a Newton iteration; for details, see Piegl and Tiller [21]. (3) Take the weighted average of n 1 and n 2 to get n. Weight by the magnitude of the parametric derivative, k dx dn k, where x is the position in physical space and kAEk denotes the Euclidean norm. That is, the new knot is given by n ¼ dx dn ð n 1 ;g 1 Þ n 1 þ dx dn ð n 2 ;g mþqþ1 Þ n 2 dx dn ð n 1 ;g 1 Þ þ dx dn ð n 2 ;g mþqþ1 Þ ðB:1Þ This weighted average produces a better splitting of an edge where the physical position changes rapidly with variation in the parameter. Most notably, this yields a better splitting of longer edges. Independent of different edge lengths, if n 1 ¼ n 2 , then n will be set to that value as well. In this way no accuracy is lost in the ''easy'' cases.  Fig. 17a. Note that uniform splitting in the parametric space translates to uniform splitting of the short edge (i.e., the circle) but quite non-uniform splitting of the longer edges opposite to it. The advocated scheme behaves in opposite fashion. The former approach produces a better mesh near the circular hole whereas the latter approach produces a somewhat more balanced mesh globally. Both approaches should be useful. The amount of time taken for the advocated refinement process is negligible in comparison with analysis. Approximate the midpoint of the element edge by P 1 , the midpoint of the chord connecting the corners of the element. Perform point inversion to find n 1 , the parameter of the point on the curve that is closest to P 1 in the Euclidean norm. Repeat the process on the opposite side of the domain to find n 2 . The knot to be inserted is a weighted average of n 1 and n 2 .