Finite Element Analysis of Structures Using C-Continuous Cubic B-Splines or Equivalent Hermite Elements

We compare contemporary practices of global approximation using cubic B-splines in conjunctionwith doublemultiplicity of inner knots (C-continuous) with older ideas of utilizing local Hermite interpolation of third degree. The study is conducted within the context of the Galerkin-Ritz formulation, which forms the background of the finite element structural analysis. Numerical results, concerning static and eigenvalue analysis of rectangular elastic structures in plane stress conditions, show that both interpolations lead to identical results, a finding that supports the view that they are mathematically equivalent.


Introduction
Structural analysis is usually performed using commercial codes that include finite elements of low (usually first or second) degree, where the accuracy of the calculations increases by mesh refinement (ℎ-version). Alternatively, keeping the number and the position of the nodal points unaltered, the numerical solution improves using polynomials of higher degree ( -version) [1].
As an extension of the above -version "philosophy, " tensor-product Lagrange polynomials as well as CAD-based (Gordon-Coons) macroelements-based on several interpolations-have been used [2][3][4]. The aforementioned macroelements integrate the solid modelling (CAD: computeraided-design) with the analysis (CAE: computer-aided-engineering). In more detail, these macroelements use the same global approximation for both the geometry and the displacement vector. In order to avoid the undesired numerical oscillations caused by Lagrange polynomials of high degree, the next generation of CAD-based macroelements replaced them with tensor-product B-splines [5]. Since 2005, the nonuniform-rational-B-splines (NURBS) interpolation has started to prevail [6].
A careful study of literature reveals that most of recent papers referring to the so-called isogeometric analysis (IGA) start with some essentials on the definition of B-splines and relevant recursive formulas due to de Boor [7]. It should be recalled that NURBS is an extension of B-splines (nonuniform and rational) modified on the basis of weighting coefficients, thus producing fully controlled sculptured surfaces [8,9]. In a B-splines expansion, the multiplicity of the inner knots plays a significant role in the continuity of the variables. In general, the multiplicity of inner knots per breakpoint in combination with a piecewise polynomial of degree ensures − -continuity of the variable (here: displacement components) [7,9]. Thus considering cubic B-splines ( = 3) in conjunction with double inner knots ( = 2), 1 -continuity is ensured. Höllig [5, page 93] has solved plane stress problems using B-splines of degree = 2, 3, 4, and 5, but his study is not a complete investigation on the influence of the multiplicity and corresponding continuity of variables involved.
On the other point of view, tensor-product Hermite elements of third degree have been proposed for the solution of fourth-order problems, such as plate-bending problems, using Galerkin-Ritz formulation. The need for smoother ( 1 ) global basis functions is also encountered in second-order problems when collocation finite element methods are utilized [10, page 66].
With these situations in mind, we next examine the relationship between particular Hermite elements of third degree 2 Journal of Structures and cubic B-splines elements with multiplicity = 2. Given that both a global interpolation (B-splines) and a local one (Hermite elements) ensure 1 -continuity, the question arises whether and how these are mutually equivalent.
Numerical examples concerning rectangular structures in plane stress conditions reveal that both approximations are equivalent to one another.

Governing Equations.
In the case of an elastic, isotropic, and homogeneous 2D structure, the governing equation in the domain Ω is given by is the stress operator, is the mass density, andü = {̈V } is the acceleration vector.
Considering the Hookean elasticity matrix in 2D problems (plane stress or plane strain), E = ( 11 12 0 21 22 the relationship between the stress, , and the strain In addition, the relationship between the strain and the dis- Therefore, the final governing equation (Naviér-Kirchhoff) is written as follows: where is a quadratic operator.

Boundary Conditions.
Let us consider the field of displacements u( , ) = { ( , ) V( , )} on the domain Ω, which is required to fulfil the governing equation (5). The boundary conditions that correspond to this problem are of two types as follows: (a) essential conditions, such as u = u( ) on Γ 1 (Dirichlet type), (b) natural conditions of the type t = t( ) on Γ 2 (Neumann type), with t denoting the traction, the components of which are related to the stress tensor as follows: where and are the components of the outward normal vector on the boundary Γ 2 . The total boundary is Γ = Γ 1 + Γ 2 .

Global and Local Interpolation
Below we present the two interpolations that will be compared to each other.

The B-Splines as a Global Functional Set.
In the evolution of time, B-splines have appeared in two different forms.
Introducing the truncated power as which is −1 -continuous, the original expression consists of a power series in the form [11,12]: It is apparent that (9) includes ( + 3) terms and ensures 2continuity, because for simplicity we considered = 3. Alternatively, (9) can be modified so as to include additional truncated polynomials of second degree; that is, Obviously, (10) includes 2( + 1) terms and ensures 1continuity.

Contemporary Procedures.
The breakthrough was made in 1972, independently by de Boor [7] and Cox [13], who both achieved the B-spline function and its derivatives to be rapidly computed through recursive formulas. The latter procedures exist today in MATLAB (spline toolbox) under the name "spcol. " In brief, for a nondecreasing set of ( + 1) breakpoints ( 0 , 1 , . . . , −1 , ) and for a certain polynomial degree " , " we define the knot vector "V": Journal of Structures 3 which highly depends on the chosen multiplicity of internal knots (usually single or double), as follows.
(i) The Cartesian coordinate is approximated in terms of control points as (ii) The variable is approximated as It is worth mentioning that the coefficients in (19) are generally different than the nodal values associated with the breakpoints, except for the ends where 0 ≡ 0 and ≡ .

Piecewise Hermite Approximation.
We refer again to the sequence 0 < 1 < ⋅ ⋅ ⋅ < , which was mentioned in Section 3.1. Let Ω = [ −1 , ] be an arbitrary element in a finite element partition of the interval [ , ] ≡ [ 0 , ]. A 1cubic element is obtained by implementing the function and derivatives at the ends −1 and of each element.
The actual element Ω is transformed linearly to the master element Ω = [0, 1] by the map On Ω, interpolating and ≡ / at the end nodes = 1, 2, the cubic Hermite expansion has the form: where the Hermite basis functions { 0 , 1 } satisfy the interpolation properties at end nodes 1 = 0, 2 = 1: for local nodal indices = 1, 2 and = 1, 2. Using properties (21) we may construct the Hermite cubics directly as

Contemporary versus Older Definitions.
For the particular case of = 3 (cubic splines) under this study, it is trivial to prove that (i) when the multiplicity equals one ( = 1), the ( + 3) basis functions involved in (19) are equivalent (not identical) to those basis functions in (9); (ii) when the multiplicity equals two ( = 2), the 2( + 1) basis functions involved in (19) are equivalent (not identical) to those basis functions in (10); (iii) therefore, a lot of research conducted in 1960s using the older framework ( (9) and (10)) may be repeated on the new unified framework of (19).
In this paper, we are concerned with = 2, which corresponds only to (10).
In order to get a better idea, the one-dimensional unit domain was divided into four subdivisions (Ω = [0, 1], = 4), and then (i) local shape functions using Hermite polynomials of third degree and (ii) global basis functions using cubic B-splines (in conjunction with = 2) are compared in and the cubic B-splines (shown in Figure 1(b)). In more detail, concerning the Hermite polynomials, those two shape functions that correspond to flexural DOFs vary between 0 and 1, whereas those two that correspond to slopes obtain both positive and negative values. In contrast, all B-splines basis functions are nonnegative and generally less than unity.

Tensor-Product Shape Functions
Let us consider a rectangular domain Ω = (ABCD) = [0, ] × [0, ] in R 2 . The axis origin is chosen at corner A, whereas the Cartesian axes and lie on the sides AB and AD, respectively. Without loss of generality, the sides (AB, CD) and (BC, DA) are uniformly divided into and segments, respectively, thus leading to ( + 1) breakpoints along AB or CD, as well as ( + 1) breakpoints along BC or DA. In this paper, we use a unique polynomial degree (cubic: = 3) for both -and -directions. As a result, the univariate function ( , 0) along the side AB can be interpolated through a piecewise B-splines polynomial of third degree in , whereas the function (0, ) along the side DA can be interpolated through a piecewise B-splines polynomial of third degree in .

Tensor-Product B-Splines.
Given the uniform subdivisions and of the intervals [0, ] and [0, ], respectively, the breakpoints along each of the four sides (AB, BC, CD, and DA) are determined. Moreover, given the multiplicity of internal knots, = 2, as well as the polynomial degrees = 3, the control points in the -and -directions are also determined. If the patch is curvilinear, then -andcoordinates have to be replaced by the -and -normalized coordinates, respectively.
While in older B-splines formulation [11,12] the degrees of freedom are associated with the ( + 1) × ( + 1) nodal points x lying at the intersections of th and th lines perpendicular to the axes and passing through the breakpoints ( , ), in this-let us say-"modern" formulation we have to deal only with the tensor-product of control points. Since the multiplicity of internal knots is = 2, the tensor-product consists of = 4( + 1)( + 1) coefficients , for each displacement component.
Therefore, the two-dimensional global shape functions are given by (the double subscript is here to emphasize the two directions) Since the multiplicity is equal to two, then the univariate approximation is 1 -continuous (whereas for 2-D ∈ 1,1 (Ω st )); Ω st = [0, 1] × [0, 1] is the standard reference square.
In general, the control points are divided into two categories, that is, ,in in the interior of the domain Ω and , (generally) near the boundary ( = ,in + , ). In more detail, if a side of the quadrilateral ABCD (e.g., AB) is straight (as in this study) the corresponding control points lie on this side (AB). In contrast, if the side is curved, then the extreme control points (P 0 and P ) will belong to the boundary and even they coincide with the corners (e.g., A and B), whereas the rest will be either inside or outside the domain Ω in accordance with the curvature of the curve AB.

Tensor-Product
where 00 = ( , ), 10 = / at ( , ), and so on. Using the maps → and → , we obtain the corresponding tensor-product shape functions on the rectangle Ω . Considering the chain rule, the derivatives in the and directions in (25) transform to partial derivatives in and .
There are now four element shape functions associated with each node to interpolate the corner values of , ≡ / , ≡ / , and ≡ 2 / . For each variable , the element has 16 degrees of freedom. The form of the element interpolant (25) transforms to where { } = { ( , ), ( , )/ , ( , )/ , 2 ( , )/ } and we have suppressed the superscript on the element shape functions for notational simplicity. Details are given in the Appendix.
This tensor-product basis ensures 1 -continuity, as is evident from the form of the global basis functions at an interface between two elements in the discretization. Equivalently, differentiating (26) with respect to and evaluating ( , ) on side = yields a cubic function ( ) on this side that is uniquely determined by the function value and its derivative ( ) at the endpoints of that side.

General.
For the given partial differential equation (5), we seek an approximate solution to (5) In this paper, the candidate functions are cubic B-splines ((23) in conjunction with generalized coefficients ) or Hermite polynomials ((25) in conjunction with kinematic degrees of freedom ). It should be clarified that the variablẽin (27) stands for either the horizontal or the vertical displacement components at any point of the elastic structure.
Based on the global shape functions involved in (27), we can apply the well-known Galerkin-Ritz method.
Without loss of generality, the boundary consists of̃1 breakpoints (which correspond to 1 control points) under Dirichlet and̃2 ones (which correspond to 2 control points) under Neumann boundary conditions. For the sake of brevity, below we limit the discussion in the two typical cases of boundary conditions, that is, Dirichlet-type and Neumanntype boundary conditions.
After imposing the boundary conditions, the dimensions of each matrix in (28) become eq × eq , where eq is the number of equations (i.e., the number of unrestricted degrees of freedom: coefficients or kinematic DOFs). For a certain choice of × subdivisions, eq has the same value either the B-splines or the Hermite formulation is implemented.

Numerical Implementation.
The numerical procedure is performed as follows. The domain is uniformly divided into a certain number of × subdivisions (cells), thus leading to ( + 1) × ( + 1) nodes.
Concerning the determination of the proper Gaussian quadrature, we start from the observation that the elements of the mass matrix are products of two basis functions, each of piecewise th (i.e., third) degree. In the particular case of a rectangular domain which is the topic of this paper, the integrant becomes of piecewise sixth degree and thus it requires four-point Gauss quadrature per direction, that is, sixteen Gauss points per integration cell.

Numerical Examples
Two MATLAB codes of very similar architecture were developed on a standard PC Pentium IV as follows: (i) a code based on rectangular 4-node (32-DOF) Hermite elements; (ii) a code based on contemporary cubic B-splines in conjunction with double inner knots; the basis functions , and their derivatives were created using the "spcol" function, which exists in the Spline Toolkit.
The eigenvalues were calculated using the standard "eig" function.
The theory is now elucidated by three examples: one example deals with static analysis while the remaining two examples deal with the eigenvalue analysis of rectangular structures in plane stress conditions. With respect to the exact solution, exact , the quality of the numerical solutioñis evaluated in terms of the relative error, which was calculated as follows: where is the moment of inertia and for a beam with rectangular cross-section and unit thickness it is given by The rectangular domain is uniformly discretized using of × subdivisions along the -and -directions, respectively, as shown in Figure 2 (for = 1, 2). For these models, Figure 3 shows that the flexural displacement at the middle of the free side BC (see Figure 2) converges to the exact solution even for a very small number of subdivisions.
Example 2 (square plate in plane stress). Consider a square plate of uniform thickness under plane stress conditions, which is fixed along its entire boundary. The following parameters were used: / = 10 4 and ] = 0.30.
In the absence of an analytical solution, a parametric analysis using ANSYS (PLANE42 elements) shows that for a uniform mesh of 90 × 90 elements the solution has practically converged. The relative errors (cf. (30)) for the first six eigenvalues are shown in Figure 4. The horizontal axis corresponds to the previously mentioned number of equations (unrestricted DOF), eq .

Discussion
Bicubic B-spline interpolation has been previously used for the finite element analysis of plates [5,15] based on Galerkin-Ritz formulation, as well as in conjunction with collocation methods for potential problems [7,10].
In the advent of the isogeometric analysis [6], the older piecewise interpolations have been substituted by knots and control points, but in some cases it may be a recast of older formulations. Within this framework, it was shown that Bsplines interpolation in conjunction with double inner knots ensures 1 -continuity. In other words, although it seems that cubic B-splines interpolation is a global approximation of both geometry and variable within the entire rectangular domain, it is equivalent of splitting the domain into × subdivisions along the -and -directions, respectively, and considering local approximation within each of the aforementioned × Hermite elements. A similar coincidence has been previously noticed in problems of one dimension, however, in conjunction with the collocation method [16].
According to Carey and Oden [10, page 68], unfortunately, there is one serious shortcoming regarding the abovementioned tensor-product Hermite elements-to retain 1continuity, the discretization must be restricted to consist of rectangles in two dimensions and rectangular prisms in three dimensions. So far, two alternatives have been proposed: (i) the use of nonconforming and (ii) simplex elements. Alternatively, it is anticipated that the tensor-product Bsplines formulation (in conjunction with double inner knots) will be applicable to problems with curved boundary, at the cost of determining the Jacobian and its inverse.
Remaining within the context of rectangular elements, higher order macroelements can be used. One possibility is 8 Journal of Structures to use similar global functions with those previously used for plate bending [17]. In addition, modified Hermite polynomials, a family of 1 -rectangular elements of increasing degree, have been proposed (see, e.g., [10, page 67]). Based on the equivalency revealed in this paper, it is anticipated that similar results will be obtained when substituting the two aforementioned alternatives with proper breakpoints and corresponding multiplicities. For example, at inner points only the displacement components can be considered, whereas at extreme points their derivatives can be considered as well. While explicit expressions have been previously derived, equivalent B-splines with proper select of multiplicities are anticipated to be applicable.

Conclusions
The Galerkin-Ritz method using contemporary tensorproduct cubic B-splines with double inner knots (global interpolation) was found to be numerically equivalent to the older 1 -continuous Hermite elements of third degree (local interpolation) when applied to rectangular domains. This finding is valid for both static and eigenvalue problems. The superiority of the contemporary tensor-product cubic Bsplines is probably that, in principle, they can work for nonrectangular structures as well. Moreover, ongoing research reveals that global B-splines interpolation can readily implement several modified Hermite alternatives such as higher degree or different number of function values and derivatives at internal and external points.

Interpolation Using Cubic Hermite Polynomials
The cubic Hermite polynomials are Renaming the corners A, B, C, and D with the local numbers 1, 2, 3, and 4, respectively, and considering the abovementioned four degrees of freedom per node for the -displacement: ( , / , / , 2 / ) , (A. 3) and another four degrees of freedom per node for the -displacement: the vector of degrees of freedom is of the dimensions 32 × 1 as follows: (A.6) and so on. Finally, for each element the interpolation of the components, and V, will be

Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.