A Review on Attempts towards CAD/CAE Integration Using Macroelements

In this paper we review several aspects of older and contemporary attempts to integrate computer-aided design (CAD: geometric model) and computer-aided engineering (CAE: finite elements, boundary elements, etc.). After a short review on formulas for the description of CAD surfaces, a systematic mechanism for creating several types of corresponding isoparametric macroelements is presented. Gordon-Coons is initially applied in conjunction with piecewise linear, Lagrange polynomials and natural B-splines. Then, it is extended to more basis and blending functions. In addition to the well-known ‘Lagrangian’-type elements, equivalent ‘Bézierian’-type elements are introduced. Tensor product B-splines and aspects of NURBS isogeometric formulation are given. In addition to quadrilaterals, triangular macroelements based on Barnhill’s interpolation are presented for the first time. The review covers applications of CAD-based macroelements in conjunction with the Galerkin-Ritz formulation, the Boundary Element Method, as well as recent Global Collocation procedures. A numerical example on a vibrating membrane elucidates the performance of the CAD-based global interpolation and depicts its superiority over the conventional finite element method.


Introduction
Computer-Aided Design (CAD) is the use of computer systems to assist in the creation, modification, analysis, or optimization of a design. One the earliest personalities who have contributed in this area, and had a vision of interactive computer graphics as a powerful design tool, were Steven Anson Coons (1912Coons ( -1979, a professor in the Mechanical Engineering Department at Massachusetts Institute of Technology (MIT) during the 1950s and 1960s. During World War II, he worked on the design of aircraft surfaces, developing the mathematics to describe generalized "surface patches". At MIT's Electronic Systems Laboratory he investigated the mathematical formulation for these patches, and in 1967 published one of the most significant contributions to the area of geometric design, a treatise which has become known as "The Little Red Book" [1]. His "Coons Patch" was a formulation that presented the notation, mathematical foundation, and intuitive interpretation of an idea that would ultimately become the foundation for other surface descriptions that are commonly used today, such as B-spline surfaces, NURB surfaces, etc. His technique for describing a surface was to construct it out of collections of adjacent patches, which had continuity constraints that would allow surfaces to have curvature which was expected by the designer. Each patch was defined by four boundary curves, and a set of "blending functions" that defined how the interior was constructed out of interpolated values of the boundaries. The interested reader may find more details in [2].
Computer-Aided Engineering (CAE) is the broad usage of computer software to aid in engineering tasks. It includes computer-aided design (CAD), computer-aided analysis (CAA), computer-integrated manufacturing (CIM), computer-aided manufacturing (CAM), material requirements planning (MRP), and computer-aided planning (CAP). In a more strict sense, CAE areas include:  Stress analysis on components and assemblies using FEA (Finite Element Analysis);  Thermal and fluid flow analysis Computational fluid dynamics (CFD);  Multibody dynamics (MBD) & Kinematics;  Analysis tools for process (manufacturing) simulation for operations such as casting, molding, and die press forming.  Optimization of the product or process.  Safety analysis of postulate loss-of-coolant accident in nuclear reactor using realistic thermal-hydraulics code.  In general, when omitting the abovementioned planning tasks, there are three phases into which any CAE procedure is divided:

62
A Review on Attempts towards CAD/CAE Integration Using Macroelements  Pre-processing -defining the model and environmental factors to be applied to it (typically a finite element model, but facet, voxel and thin sheet methods are also used).  Analysis solver (usually performed on high powered computers).  Post-processing of results (using visualization tools). This cycle is iterated, often many times, either manually or using built-in procedures, or finally by linking the entire procedure into a commercial optimization software such as [3,4].
Henceforth, CAE will restrict to the finite element (FEM/FEA) and the boundary element (BEM) analysis, whereas a very recent Global Collocation Method (GCM) will be reviewed and commented. It is reminded that FEM decomposes the model of the structure in small areas or volumes, called finite elements [5,6], whereas BEM deals mostly with the boundary of the structure [7].
The decade of 1960s was very important for designers, as since then they have long used computers for their calculations. Initial developments were carried out in the 1960s within the aircraft and automotive industries in the area of 3D surface construction and NC programming, most of it independent of one another and often not publicly published until much later. Some of the mathematical description work on curves was developed in the early 1940s by Isaac Jacob Schoenberg, Apalatequi (Douglas Aircraft) and Roy Liming (North American Aircraft), however probably the most important work on polynomial curves and sculptured surface was done by Pierre Bézier (Renault), Paul de Casteljau (Citroen), Steven Anson Coons (MIT, Ford), James Ferguson (Boeing), Carl de Boor (GM), Birkhoff (GM) and Garabedian (GM) in the 1960s and W. Gordon (GM).
It is argued that a turning point was the development of SKETCHPAD system in MIT in 1963 by Ivan Sutherland [13], who was a student of S.A. Coons; however his PhD thesis was supervised by C.E. Shannon (1916 -2001), the "father of information theory". The distinctive feature of SKETCHPAD was that it allowed the designer to interact with his computer graphically: the design can be fed into the computer by drawing on a CRT monitor with a light pen. Effectively, it was a prototype of graphical user interface (GUI), an indispensable feature of modern CAD.
First commercial applications of CAD were in large companies in the automotive and aerospace industries, as well as in electronics. Only large corporations could afford the computers capable of performing the calculations. Notable company projects were at GM (Dr. Patrick J.Hanratty) with DAC-1 (Design Augmented by Computer) 1964; Lockhead projects; Bell GRAPHIC 1 and at Renault (Bézier) UNISURF 1971 car body design and tooling. One of the most influential events in the development of CAD was the founding of MCS (Manufacturing and Consulting Services Inc.) in 1971 by Dr. P. J. Hanratty [14], who wrote the system ADAM (Automated Drafting And Machining) but more importantly supplied code to companies such as McDonnell Douglas (Unigraphics), Computer vision (CADDS), Calma, Gerber, Autotrol and Control Data. As computers became more affordable, the application areas have gradually expanded. The development of CAD software for personal desk-top computers was the impetus for almost universal application in all areas of construction. Other key points in the 1960s and 1970s would be the foundation of CAD systems United Computing, Intergraph, IBM, Intergraph IGDS in 1974 (which led to Bentley MicroStation in 1984) CAD implementations have evolved dramatically since then. Initially, with 2D in the 1970s, it was typically limited to producing drawings similar to hand-drafted drawings. Advances in programming and computer hardware, notably solid modeling in the 1980s, have allowed more versatile applications of computers in design activities.
Key products for 1981 were the solid modelling packages -Romulus (ShapeData) and Uni-Solid (Unigraphics) based on PADL-2 and the release of the surface modeler CATIA (Dassault Systemes). Autodesk was founded 1982 by John Walker, which led to the 2D system AutoCAD. The next milestone was the release of Pro/ENGINEER in 1988, which heralded greater usage of feature-based modeling methods and parametric linking of the parameters of features. Also of importance to the development of CAD was the development of the B-rep solid modeling kernels (engines for manipulating geometrically and topologically consistent 3D objects) Parasolid (ShapeData) and ACIS (Spatial Technology Inc.) at the end of the 1980s and beginning of the 1990s, both inspired by the work of Ian Braid. This led to the release of mid-range packages such as SolidWorks in 1995, SolidEdge (Intergraph) in 1996, and IronCAD in 1998. Nowadays CAD is one of the main tools used in designing products.
In order to analyze mechanical components, which are usually characterized by complex shapes, the computer-aided-design (CAD) model is usually exported in the form of a neutral file (IGES, DXF, STEP) that is further imported by a selective finite element (FEM) or boundary element (BEM) code. This procedure induces some difficulties such as the appearance of double points (or double lines) or loss of at least a few geometrical data. Of course, specific CAD converters such as PATRAN can export many types of FEM data that can be immediately feed commercial FEM codes. Specialized CAD/CAM/CAE integrated systems such as PRO/ENGINEER, IDEAS and Dassault Systèmes (SOLIDWORKS, PLM: Product Lifecycle Solutions) exist for the last 20 or 25years. However, in all these tools the computational mesh generation is a time-consuming task that has to follow the creation of the CAD model. A survey of thirty-four finite element systems until 1981 was reported by Brebbia [15].
The term "CAD/CAE integration" may have multiple meanings, which will be better understood through an example. Suppose that a bicycle designer assigns a tentative tube size for the skeleton of the structure (CAD); then, an integrated CAD/CAE system has to automatically understand not only the diameter and the thickness but also the second moment of inertia of this tube, without the engineer having to perform the relevant numerical calculations. Closely related, the engineer can generate the mesh without exiting the CAD system and then entering the CAE system; in other words, both systems are integrated in one, sharing the same database as is possible.
A second point of view is the usually reported fact that the CAD/CAE integration and scientific visualization is related to both the reduction of errors in the information transfer from system to system as well as the reduction of memory resources and overall computational time [16,17]. The 'vehicle' for CAD/CAE integration is to use a common platform to describe the geometry and the unknown variable that satisfies a partial differential equation (PDE) of a multi-physics boundary-value-problem (e.g. heat transfer, elastostatics, thermo-elasticity, fluid mechanics, electromagnetic, etc). A third point of CAD/CAE integration is to replace the usual finite elements (of small-size) with others of larger size applying the same interpolation for both the geometry and the variable. Relevant macro-elements are usually of "isoparametric" type whereas "isogeometric" ones have appeared [18].
The lack of a complete, synthetic and critical report on the abovementioned third point of CAD/CAE integration, which is directly related to the development and use of macroelements, has motivated the author to write this review article.
The paper is structured as follows. Section 2 summarizes the one-dimensional (univariate) interpolation as well as the five basic CAD formulations. Section 3 reviews the usual applications in which Coons' interpolation has been used, as well as the derivation of global shape functions. Section 4 provides details about the isoparametric or isogeometric-like approximation that is used for the interpolation of both the geometry and the variable or the unknown coefficients. Section 5 reviews the three possible computational methods to be applied in conjunction with the global interpolation within subregions, and provides the numerical procedure for the estimation of mass and stiffness matrices of the structure. Section 6 presents original numerical results concerning the extraction of natural frequencies for a fixed rectangular membrane. Section 7and Section 8are the "Discussion" and the "Conclusions", respectively.

General
Historically, there are rather five main interpolation formulations in CAD theory [19]: (i) Coons interpolation, (ii) Gordon-Coons interpolation, (iii) Bézier interpolation, (iv) B-splines interpolation, and (v) NURBS. The first two of them directly refer to the coordinates of nodal points that belong to the boundary or the interior of a patch, whereas the last three refer to control points (at the edges of polygonal generators that generally do not belong to the curvilinear boundary). In all five cases the interpolation of geometry refers to the coordinates of either nodal points or control points multiplied by proper basis functions. For the sake of completeness, we start with one-dimensional and then extend to two-and three-dimensional interpolation.

One-Dimensional Interpolation
Let us consider the domain a≤x≤b in which we wish to interpolate the function f(x).

Power Series
The function is approximated by:

Lagrange Polynomials
The function is approximated by: denotes the well-known Lagrange polynomial (see, for example, [20]).

P-Method
The relevant series expansion includes the nodal values at the ends, f(a) and f(b) for which the classical linear hat functions are considered, i.e.: It is reminded that the Legendre polynomials are given by and they are solutions of the ordinary differential operation: It is also reminded that the well-known Gauss points are roots of the above Legendre polynomials.
The advantage of these functions is that they are orthogonal thus leading to banded mass and stiffness matrices [21, pp. 37-47].

2.2.4.ChebyshevPolynomials
Although these polynomials do not directly appear in CAD formulations, it is instructive to mention that they can approximate a function with the same accuracy with the power series. They are categorized in polynomials of first and second kind. For example, the Chebyshev polynomials of the second kind are defined as: and their roots are given by As we will see later, the above roots can be used as collocation points in a global collocation method, which is an alternative to the Galerkin-Ritz procedure [22,23].
The geometrical coefficients, { } i P , are called control points. Equation (10) denotes that for a curve that is determined by (n + 1) control points, the term of highest degree is n u . This implies that the polynomial degree of Bézier curve is determined by the number of control points. The basic properties of a Bézier curve are many and can be found elsewhere [24][25][26][27]. In brief:  The curve passes through the first (Ρ 0 ) and the last (Ρ n ) control points and its tangent at these ends has the direction of the first straight segment (Ρ 0 Ρ 1 ) and the last one (Ρ n-1 Ρ n ), respectively, of the control polygon (Ρ 0 Ρ 1 … Ρ n-1 Ρ n ).  They are rapidly calculated by the following recursive formula:  Their derivative is also calculated by the following recursive formula: it became necessary to introduce the n-th rational Bézier curve, which is defined as: The meaning of splines was published in 1946 and later by Schoenberg [28]. It refers to the points n n x f x which we wish to interpolate through a multiply-defined function ( ) f x . The points 0 1 , , , n x x x  are called "breakpoints". For the sake of briefness and without loss of generality, we reduce to (cubic) polynomials of third degree. The desired properties are as follows:  In each interval the function ( ) f x is a cubic polynomial.
 The function ( ) f x , as well as the first and second derivatives, are continuous at the above points. Introducing the truncated power as: which has C (m-1) -continuity, the original expression consists of a power series in the form [28]: ( ) It is apparent that (18) includes (n+3) terms and ensures C 2 -continuity, because for simplicity we considered m = 3. Alternatively, (18) can be modified so as to include additional truncated polynomials of second degree, i.e.: ( ) Obviously, (19) includes 2(n+1) terms and ensures C 1 -continuity.

Contemporary Procedures
The breakthrough was made in 1972, independently by Cox [29] and DeBoor [30], who both achieved the B-splines and their derivatives to be rapidly computed through recursive formulas [including (11) and (12)]. Charles DeBoor [26] continued his research and his initial FORTRAN codes exist today in MATLAB (spline toolbox) under the name 'spcol', among others [31].
In brief, for a nondecreasing set of (n+1) breakpoints , and for a concrete polynomial degree 'p', we define the knot vector 'U', which strongly depends on the multiplicity of the internal nodes, as follows: Based on either of the knot vectors of (20), (m+1) control points (Ρ 0 Ρ 1 … Ρ m-1 Ρ m ) are constructed, so as the following relationship always holds:

Contemporary Versus Older Definitions
Some young readers may believe that (22) is different than (18) or (19). However, for example in case of p = 3 (cubic splines), it is trivial to prove that:  When the multiplicity equals to one, the (n+3) basis functions of (22) are equivalent (not identical) to those basis functions given by (18).  When the multiplicity equals to two, the 2(n+1) basis functions of (22) are equivalent (not identical) to those basis functions given by (19).  When there are no internal points (n = 1), B-splines coincides with p-thBézier.
 Therefore, a lot of research conducted in 1960s using the older framework ( (18) and (19)) may be repeated on the new framework of (22).

2.2.7.NonUniform Rational B-Spline (NURBS) Curve
Similarly to the rational Bézier curve [see, (16)], non uniform rational B-splines (NURBS) is an extension of the B-splines in the following form: where {w i } are the weights. Usually, a = 0, b = 1 and w i > 0 for all i. Taking the rational basis functions as, the curve is expressed as: (25) and they are piecewise rational functions in the interval

Comparison between Different Sets Of Basis Functions
It is well known that power series is identical with Lagrange polynomials [24, pp. 304-308]; the only difference is that the first includes arbitrary coefficients, whereas the second cardinal shape functions (of [1,0]-type) that are multiplied by the associated nodal values.
Chebyshev polynomials are also multiplied by arbitrary coefficients and they may have some advantages when collocating at their roots [22,23].
The use of Legendre polynomials is somehow different. The relevant series expansion includes the nodal values at the ends, f(a) and f(b) for which the classical linear hat functions are considered, i.e.: In addition, a certain desired number of "bubble" functions, in the form of differences between Legendre polynomials, are usually used. The advantage of these functions is that they are orthogonal thus leading to banded matrices [21, pp. 39-46]. Bernstein polynomials have been used in the construction of the well-knownBézier curves, and can be found in many textbooks such as [25][26][27].
B-splines is nothing new but the usual spline interpolation, which can be found in every textbook of numerical analysis such as [20], whereas detailed information is given elsewhere [27][28][29][30]. The superiority of B-splines is that it can handle every high number of data, in contrast to Lagrange polynomials that cannot handle more than ten to twelve points at maximum.
Finally, NURBS is an extension of B-splines as it uses weights for every control point.
As one can notice in Fig.1, while Lagrange polynomials may exceed the unity and be characterized by high values particularly in the neighbourhood of the ends, Bernstein polynomials (or even general B-splines, not shown) are always less than or equal to unity (at the end nodes). It is worth-mentioning that, in both cases, the sum of the shape functions is equal to unity (rigid-body property, partition of unity). Similarly, in the p-method the sum of the first and last basis functions (N 1 and N 2 , associated to the ends of the domain) is equal to unity, whereas the intermediate "bubble functions", j φ , vanish at the ends but their sum is different than zero at every internal point (see, [21]). A side ascertainment on the equivalency between the three functional sets in power series, i.e. Lagrange polynomials, Bernstein polynomials and the p-method is the numerical finding that the calculated eigenvalues are identical [32,33].

Two-Dimensional (2D) Interpolation
In brief, surfaces are categorized into quadrilateral and triangular patches, i.e. with four and three sides, respectively. The interpolations can treat either straight or curvilinear sides of the patch. The oldest way to create a 2D patch is Coons interpolation [1]. Later, tensor products of Bézier, B-splines and NURBS were applied [19,27].

Quadrilateral Patches
We distinguish two categories. The first deals with data along the four boundaries of the quadrilateral; in other words, it consists the "boundary-only" formulation with nodal points along the boundary only. The second formulation uses the aforementioned boundary data as well as additional data in the interior of the quadrilateral; therefore it includes internal nodes as well, and is called "transfinite" interpolation. can easily be mapped to a unit square in the ξη parametric domain shown on the right of the same figure by the method of Coons' patch [1]. For purposes of generalization, the relevant theory is given below using suitable projections. First, the concept of the lofting projector 'P' is introduced. This projector is any idempotent linear operator, which maps a true surface to an approximate surface, subject to certain interpolatory constraints.
Let us assume that the Cartesian co-ordinates x in A, with ξ and η denoting normalized co-ordinates, are known at the boundaries ( 0,1 ξ = ; 0,1 η = ) of a curvilinear patch of area A. Let us also define the well-known cardinal blending functions: The above lofting operators form the basis for the definition of more complex operators with blending interpolation properties in more than one direction. So, the two-dimensional lofting operator: (28) can be constructed with the aid of the unidirectional operators { } P ξ x and { } P η x . Finally, the co-ordinates of any point in the interior of the curvilinear patch is approximated as: or, using conventional notation, as: According to Gordon's interpolation formula [34], the Cartesian coordinates of each point P of in the patch can be approximated by its as well as by its inter-boundaries (if any) as follows: For example, based on Gordon's interpolation, it is easy to derive the shape functions for a five-node element having one internal noded at its centre (in addition to the four ones at the corners) [35,36]. Remark: It has been found that when the nodal points on the opposite sides of a patch possess identical normalized positions with those nodes in the interior, then Gordon-Coons interpolation degenerates to the cross product P ξη [37, p.328]. Under the latter conditions, when Lagrange polynomials have been used, Coons interpolation coincides with the well-known interpolation in finite elements of Lagrange type (or Lagrange family) [5, p.153]. In conclusion, in patches where structured meshes of nodal points are used, all three basic CAD interpolations (Coons, Bézier, and B-splines) are expressed by a tensor product of 1-D ξ-and η-interpolations. Details are given below. By defining the breakpoints along each side of a curvilinear quadrilateral ABCD or a box-like volume ABCDEFGH, for a certain polynomial degree (p) and for a certain multiplicity (usually one or two), the basis functions ( ) , i p N u are determined. It is noted that here u is either of the normalized coordinates ξ, η, or ζ, whereas the number of basis functions equal to the (m+1) control points in the particular direction (m ξ , m η , and m ζ , in the ξ-, η-, or ζ-direction, respectively). Then, the aforementioned one-dimensional control points produce n = (1 + m ξ )×(1 + m η ) or n = (1 + m ξ )×(1 + m η )×(1 + m ζ ) control points for 2D and 3D problems, respectively. Each of the aforementioned control point is associated with one global shape function given by: Therefore, the geometry is interpolated by the expressions:

Triangular Patches
We will present two alternatives as follows.

Side Degeneration
One possibility is to use the formulas previously used in quadrilateral patches, by degenerating one of the four sides into one point [1, p.3]. For example, when degenerating the entire side AD into a unique point A, the quadrilateral degenerates into a triangle ABC, (30) The mathematical description of triangular patches has been published by Barnhill [39][40][41]. In those works a "trilinearly blended interpolant" formulation is introduced, which is quite different than Gordon-Coons; today this formula can be found even in CAD textbooks such as [42, p.244]. The triangular patch shown in Fig. 4a can be mapped into the parametric patch of Fig. 4b through the relationship: The parametric patch of (36) is expressed as: In this case, the parametric patch can be intersected with increasing values of the variables u,v that vary between 0 and 1, and then calculating the corresponding values of w from each set of u and v.
Apart from CAD applications, no application of Barnhill's formula to engineering analysis has been reported so far. For the first time, it is shown that, (36) inherently includes the conventional three-and six-nodedtriangular finite elements. Details are given in Appendix A and Appendix B, respectively. Of course, (36) allows for the derivation of novel arbitrary-noded triangular macroelements, that have been successfully applied to stress analysis problems [90].

Three-Dimensional (3D) Interpolation
Volume blocks are categorized into hexahedral (eight vertices, and twelve edges) and tetrahedral (four vertices, and six edges) blocks. The interpolations can treat either straight or curvilinear edges of the volume.

General
A boxlike region ABCDEFGH, shown in Fig. 5a, can easily be mapped to a unit cube in the rts parametric domain (0 ≤r,s,t≤ 1) shown in Fig. 5b. The relevant formula may be found in [1, p.41] and, as previously mentioned, is has been applied by Cook [43] for finite element mesh generation purposes. For more details the interested reader may also consult a textbook [44]. The only difference with a quadrilateral patch is that, here, six functions (associated to the six faces) are blended, instead of four functions (associated to the four boundary edges of the 2D patch). Below, the same formula is written in terms of projections.
Besides the above-mentioned Again, summation over repeated indices is understood.
Having introduced the one-, two-and three-dimensional operators, then the following formula describes the interpolation of the co-ordinate vector [1,p.41]:   (39) is generally applicable but it can be further simplified and be written in more readable expressions. So, in the case of a generalized curvilinear hexahedral /paralleloid (boxlike region), the geometry includes:  Six surfaces, S  Twelve edges, E, and  Eight corners, C Obviously, (39)includes all three quantities: Surfaces (S), edges (E) and corners (C). For convenience, the projections related to the S, E and C are denoted as follows: However, in the case of adequately smooth and regular surfaces, the edges E can sufficiently describe S. In fact, applying (31) on the six surfaces of the superbrick, one can easily derive the following relationship: Then, substituting (40) and (41) into (39), one can finally derive three equivalent expressions of the three-dimensional Coons' interpolation formula, as follows:

Comments on 3D Equivalent Expressions
Obviously, (42) is the most general expression because it includes any type of the surrounding surfaces S (edge-only noded and/or internally noded). Moreover, (43) is obtained by eliminating the edges (E) and it is based on the surface S, the last being corrected by the co-ordinates of the corners C. Finally, (44)includes only the twelve edges E and eight corners C, or in other words, the absolutely necessary data for the construction of a Coons' block made of Coons' surfaces; it is evident that (44) uses nodal points arranged along the twelve edges only. In conclusion, in cases where the superbrick is sufficiently regular, (44) is the most advantageous. Using conventional notation, (44) becomes:

The Three Alternative Uses of Coons Interpolation
In more detail, bivariate (2D) Coons interpolation on curvilinear surfaces is useful for (i) mesh generation for two different purposes (2D FEM or 3D BEM applications), and (ii) to derive 2D finite macroelements.
It has been previously reported that, the mesh produced by 2D-Coons interpolation can be easily smoothed by taking the average of the eight neighbouring nodes. Also, the mesh produced by 3D-Coons interpolation can be easily smoothened by taking the average of the twenty-six neighbouring nodes [45][46][47].

Older Findings
From the analysis sited in Section 2, it becomes clear that in 1967 Coons' interpolation was the first achievement in CAD surface theory, primarily aiming at describing surfaces or volumes [1]. In early 1970s, FEM analysts understood that the same interpolation might have a second usage, as a generator of structured finite element meshes [43,[48][49][50][51][52][53][54]. In more details, first a mesh generation is constructed on the parent unit square or cube for 2D and D problems, respectively. Second, nodal points are constructed along the boundaries of the virtual structure. Third, (30), or (39) and its variations, are applied and thus establish a mapping between the aforementioned nodes in the parent element and the virtual structure. As a result, the desired finite element mesh is generated.
Moreover, it was early understood that Coons interpolation could be used for engineering analysis purposes as well [55,56].

Contribution on Conventional Finite Element Meshes
In subsection 2.4.3 we mentioned that, due to the fact that the mesh generated using (30) or (39) for 2D and 3D problems, respectively, is usually of low quality, smoothing schemes have been successfully applied [57,[45][46][47]. The convergence leads to almost rectangular elements.

Construction of Gordon-Coons Macroelements for
Engineering Analysis

Literature Survey
Early thoughts on the construction of large elements were made by the pioneer Bruce Irons [58]. Expressions of lower order large elements based on Coons-Gordon ideas are [55,[59][60][61]. Numerical examples on 2D potential problems were presented in 1989 by Kanarachos and Deriziotis [62] at the National Technical University of Athens (NTUA); as a member of NTUA's FEM/BEM group, the author can confirm that relevant attempts started therein in 1982. This work was extended to 2D elasticity [63]; in the relevant conference the author announced that the same methodology is not applicable only to Coons but also in conjunction with all other CAD formulations such as Gordon-Coons, Bézier, B-splines and NURBS.
Boundary-only Coons interpolation was applied to 2D potential problems [37,[64][65][66][67][68], axisymmetric elastostatics [69], and 2D elastodynamics [70][71][72][73]. In the process of the research it was found that boundary-only Coons interpolation corresponds to a series expansion of the solution with those terms appearing in Pascal's triangle with a surplus of two (Serendipity family), a fact that was later found as a note in [74].
In some examples, the boundary-only Coons formulation was more accurate not only than the conventional FEM, with the same mesh density, but also than BEM [75].
It was later understood that, it is not always possible to accurately solve all engineering problems unless a sufficient number of internal nodes are used. The latter is easily achieved using Gordon-Coons (transfinite) interpolation [35-37, 72, 76]. Boundary-only Coons interpolation was successfully applied to 3D potential and elasticity problems [77][78][79]. The theme definitely closes with [80], which compares twelve alternative models one another; although it focuses on 2D elastostatics, the conclusions are of general importance.
The overall conclusion is that boundary-only formulation sometimes cannot converge to the accurate solution, whereas the Gordon-Coons [88] always achieves it [89].
Concerning structures of triangular shape, the work of Kanarachos and Dimitriou [38,89] has shown reliable results when a side of the quadrilateral is degenerated. In addition to the latter works, Provatidis and Antoniou [90] have found encouraging results when Barnhill's interpolation [39][40][41] is used for the engineering analysis.
As one can see inspecting [91][92][93][94][95][96][97][98][99][100][101][102][103][104][105][106][107][108][109], Pierre Bézier was an employee of Renault (in France) and started publishing since 1966, which is almost the same period in which Steven Coons was working at MIT (in USA). Bézier's theory is related to the discovery of the "control points", which are vertices of a polygon (generator) that controls a curve. The same theory is also extended to curved surfaces, using a set of generators [19,27]. Extension to volumes is possible [110]. In the context of engineering analysis, it has been found that Bézier formulation is equivalent to the use of monomials or Lagrange polynomials [111], a matter that has been also treated in 1D problems [32,33,112].
Later, rational Bézier curves [113][114][115], which have several advantages over polynomial Bézier curves, appeared. Clearly, rational Bézier curves provide more control over the shape of a curve than does a polynomial Bézier curve. In addition, a perspective drawing of a 3D Bézier curve (polynomial or rational) is a rational Bézier curve, not a polynomial Bézier curve. Also, rational Bézier curves are needed to exactly express all conic sections. A degree two polynomial Bézier curve can only represent a parabola, in a similar way with the Lagrangian polynomials. Exact representation of circles requires rational degree two Bézier curves. Rational Bézier curves permit some additional shape control: by changing the weights, you change the shape of the curve. It is possible to reparametrize the curve by simply changing the weights in a specific manner.
Finally, in the process of CAD evolution, the use of B-Splines was proposed in 1973-1974 [116,117]. It is reminded that B-splines is nothing different than usual splines but they are replaced by non-cardinal basis functions, ( ) . i p N ξ of polynomial degree p, which are calculated for every normalized coordinate, ξ , through fast iterative formulas developed independently by Cox [29] and DeBoor [30] in 1972, as mentioned by the latter researcher in his monograph [26]. B-splines appear a local support and include, as a special case, the abovementioned Bézier curves. Despite the local support of B-splines, a better control on regional changes is achieved when using the non-uniform rational B-splines (NURBS) that were developed in mid-seventies till early eighties [118,119]. The interested reader may also consult a later survey [120] and a useful monograph of 646 pages [27].
According to the abovementioned sequence in CAD evolution, i.e. Gordon-Coons Bézier B-splines NURBS, a relevant contribution in structural (generally in engineering) analysis is anticipated, as follows. The common characteristic of all relevant works in engineering analysis is that, they all aim at decomposing the domain into CAD-controlled subregions [121][122][123][124].
In this paper we will use the term "isogeometric" when referring to interpolations related to "control points", i.e. Bézier, B-splines and NURBS.
Concerning Bézier representation, it is reminded that it is equivalent with the use of Lagrange polynomials [111]. In other words, for any selected degree of the polynomials involved, despite the fact that the control points sweep the entire domain, in fact it is "as if nodal points along only the sides"were used (2D: 4 sides of the quadrilateral, 3D: 12 edges of the hexahedron).
Concerning the evolution in NURBS, the interested reader may consult [142], among others. One of the first attempts to incorporate NURBS in structural analysis, and particularly in shape optimization, is due to Schramm and Pilkey [143], in 1993, followed by other investigators such as [144][145][146]. In 2005, the group of Prof. T. J. R. Hughes introduced the term "isogeometric" elements [147][148][149], and later they applied an efficient quadrature to reduce the CPU-time [150].
The author [68] had timely discovered that, despite its elegance, the global interpolation is generally slower than the common finite element method (Galerkn-Ritz), of course for the same mesh density. In 2005, he proposed to preserve the CAD-based global interpolation and replace the domain integration through a global collocation scheme [73, p. 6704].
Concerning the Global Collocation Method (GCM), an initial contribution is due to Blerk and Botha [151]. The problem was systematically tackled starting from one-dimensional problems, first starting with Lagrange polynomials [22,23,152] and then continuing with B-splines [153,154]. Then, the same idea was extended to 2D potential problems [155,156], elastostatics [157], as well as eigenvalue problems in acoustics and elastodynamics [158,159]. It is worth-mentioning that the application of isogeometric collocation methods was later proposed by others [160].

Global Shape Functions
Following [80], let us assume that the sides AB, BC, CD and DA of the patch ABCD shown in Fig. 3 include q 1 , q 2 , q 3 and q 4 nodes, respectively. In this case, the corresponding number of subdivisions per side are n 1 =q 1 -1, n 2 =q 2 -1, n 3 =q 3 -1 and n 4 =q 4 -1, respectively. Then, the number of nodes along the boundary of the patch becomes: q q n = + (47) In the sequence, the boundary values of the displacement vector, i.e. u(ξ,0), u(1,η), u(ξ,1) and u(0,η) in (30)   Similar interpolations are also considered along the inter-boundaries (if any). Finally, if all these interpolations are substituted into (48), and subsequently into (30) or (32) with j u denoting the displacement vector at nodal point 'i', appearing at the boundaries or/and the interior of the Coons macro-element. In case of Bézier B-splines or NURBS, j u in (49) should be replaced by j a . For example, the global shape function associated with the corner node A is given by: Also, the global shape function associated with an intermediate node along the side AB is given by: Finally, the global shape function associated with an internal node is given by:

The Role of Blending Functions
In addition to the index "i" that determines the basis functions along each of the four sides of a patch ABCD, now we deal with the type of the blending functions.
When the blending functions are linear functions influencing the entire side: we assign them the kodej = 0. This case corresponds to the "boundary-only" Coons macroelement, under the notation "Ci0". Again, "Ci0" means that all four sides are interpolated with the same set of basis functions, i.e. i = 1 for piecewise linear, i = 2 for B-splines, i = 3 for Lagrange polynomials, i = 4 for Bernstein (Bézier) polynomials, and i = 5 for NURBS. In addition to the above linear functions with compact support between two adjacents nodes (kodej = 1), (iii) cardinal cubic B-splines (kodej = 2), and (iv) Lagrange polynomials (kodej = 3). Based on the above definitions and code pairs (i, j), we can construct 5×4 = 20 macroelement models, denoted by "Cij", i = 1, 2, 3, 4, and 5, whereas j = 0, 1, 2, 3, and 4 (C stands for "Coons"). Again, the first index, i, refers to the type of the basis functions, whereas the second index, j, refers to the type of blending functions.
It is easy to understand that the number of possible models Cij can be highly increased when some

Isoparametric and "Isogeometric-like" Approximations
Gordon-Coons and Lagrangian type elements are based on the following series expansion: which characterizes all "isoparametric" elements. In contrast, "isogeometric-like" approximations (i.e., tensor products of Bézier, B-splines and NURBS) are based on the following series expansion: where j a are not generally the nodal values at the control points but generalized coefficients. Only at the very ends of the boundary the extreme control points (P 0 and P n ) coincide with the corresponding ends of generator curves and, therefore, the coefficients 0 a and n a coincide with the values 0 U and n U , respectively.
It is remarkable that in all "isogeometric-like" approximations the key problem is to impose the boundary conditions. This difficulty is not apparent in 1D problem, but it appears in 2D and 3D ones. For example, if the entire side AB of the quadrilateral ABCD is given homogeneous Dirichlet conditions (U = 0), this implies that all coefficients j a in (54) must vanish. In contrast, when half of the side AB is given the zero value whereas the rest is given for example by U = ξ-1/2, we cannot 'a-priori' calculated these coefficients, unless we fulfill the condition at so many points as the number of control points along AB, and then invert the produced matrix. In this context, according to Natekar et al. [124], bivariate NURBS representation is applied to derive shape functions ( ) η ξ , I N that are based on the set of I-th control points defining the system geometry. The same shape functions are also used to approximate the displacement field within the domain. It is remarkable that at any point within the domain the sum of these shape functions equals to the unity but the value of ( ) η ξ , I N at the I-th control point is not unity to that node. Therefore, "even if the control point were to be coincident with the location of the boundary condition, direct application of the boundary condition is not possible since the specified field value will be distributed to control points influencing the point under consideration" [124].
The abovementioned difficulties were the reason that the author preferred to work mainly with Coons interpolation (although very early he suggested the application of NURBS [45, p.1 and pp. [32][33], while other research teams directly followed the tendency of the modern CAD systems and worked with NURBS [18], among others.

General
The vibration of an elastic solid structure is described through the stress equilibrium equation [5]: while for acoustical cavities (or membranes) it holds:

Numerical Procedure
In this work we consider three numerical procedures: (i) Galerkin/Ritz, (ii) Collocation, and (iii) Bounday Element Method, as follows.

Galerkin/Ritz
Applying the Galerkin method to (55) or (56), one obtains the well-known matrix formulation [5]:  (58) where i N denote the shape functions by which the solution (displacement component or acoustic pressure) is interpolated as follows: where e n is the number of DOF per element.

Collocation
Using the same global shape function as above (section 4.1), and fulfilling the PDE (55) or (56), one obtains a similar matrix formulation with (57) with the difference that now no domain integral appears, that is: in elasticity , 1 , in acoustics , in elasticity , in acoustics where 'i' corresponds to the collocation points and 'j' to the control points.

Boundary Element Method (BEM)
In the case of elasticity problems, the co-ordinate vector within the l-th patch, possessing q l nodes, is interpolated on the basis of the boundaries of the patch as follows: Since each patch is considered as an isoparametric element it By substituting (62) in the usual integral equation (see for example [7]) and summarizing over the N p patches in which the boundary is divided, one obtains: In (63) the infinitesimal area dΓ is given by: where the Jacobian is calculated as: Now, for the purposes of the numerical integration only, the patch is divided into s r N N × cells where a second set of normalized co-ordinates ( ) s r is introduced [84]. Therefore, the term , which requires a trivial (e.g., 2 2 × , 3 3 × , 4 4 × ) Gaussian quadrature. A selective integration scheme has been recently developed.
Therefore, the final algebraic system obtains the form: (66) where U is the displacement vector of all nodes on the boundary of the structure (along the patch edges), ip U and ip P are displacement and traction vectors referring to the ip-th patch as well as ip H and ip G are the nonsymmetric influence matrices.
By properly assembling the submatrices in (66) we obtain: If q is the total number of nodes along all edges of the patches on the boundary of the structure, the dimensions of the vectors and matrices in (67) are as follows:  The above symbol q is larger than q ( q q q ∆ + = ) with Δq depending on the number of sharp corners and their multiplicity in traction discontinuity [7].
The above static analysis described by (67) can be also extended to the solution of dynamic problems too. Briefly, using a set of radial basis functions a mass matrix is constructed for the entire structure and it is combined with the static matrices H and G shown in (67) [86,87].

Numerical Results
In addition to previous results, in this section we present additional results that elucidate the computational performance of the proposed techniques.
Example: Eigenvalues of a fixed rectangular membrane The dimensions of the membrane in the x-and y-directions are a×b = 2.5×1.1 m, the wave velocity is taken as c = 1m/s, whereas the entire boundary is considered to be fixed. The free vibrations are governed by (56) and the exact eigenvalues are given by: The rectangular domain,2.5×1.1, is uniformly divided into n x ×n y subdivisions, preserving the same ratio, i.e. n x ×n y = 4×2, 6×3, 8×4, 10×5, 12×6, 14×7, and 16×8, thus leading to (n x +1)(n y +1) = 15, 28, 45, 66, 91, 120, and 153 nodal points, respectively. In all methods the number of unknowns is the same, whereas 2(n x +n y ) = 12, 18, 24, 30, 36, 42, and 48 nodal points belong to the boundary and have to be eliminated. Thus, in the particular case of Dirichlet problem, the number of unknowns is common for all the three methods as follows:  Conventional finite element (FEM) using four-node bilinear elements.  (Gordon)-Coons Patch Macroelement (CPM), in conjunction with the Galerkin-Ritz formulation (domain integration), according to (58).  The abovementioned Gordon-Coons approximation, in conjunction with the Global Collocation Method (GCM), that is without domain integration, according to (60). The computational error was determined by the formula

Discussion
It was shown that all five CAD interpolations, i.e. Coons, Gordon, Bézier, B-splines and NURBS are capable of producing global shape functions that influence the entire patch of volume, thus constructing reliable macroelements that can deal with a large portion (subregion) of a structure.
The great amount of relevant work performed by the FEM/BEM group at NTUA started in 1982 but it was published for the first time in 1989 [62]. It was initially based on the nodal point-based variable u, using natural cardinal (of [1,0]-type) cubic B-splines of Schoenberg's style without compact support (see, for example, [66] and a graph in [69, p.541]). This choice was rather in contrast to the mainstream 76 A Review on Attempts towards CAD/CAE Integration Using Macroelements investigations by other researchers, who used to deal with generalized coefficients, a . Here is the proper point to clarify that, the motivation of using Schoenberg's formula is the fact that the latter deals directly with the nodal values uthemselves, and not with arbitrary coefficients α related to the control points. However, in the particular case of cubic interpolation, the natural B-splines (used in [62][63][64][65][66][67][68]) lack the cubic monomial 3 x , a fact that can be technically recovered when using the two rotations at the ends of every of the four boundary sides [80, p. 958].
Later, it was decided to leave cardinal B-splines and work in conjunction with Lagrange polynomial interpolation along each of the four sides of the Coons' patch [89,37,72,73], among others. Lagrange polynomials have the advantage that they are directly associated to the nodal values. In this case, the results were extremely accurate and the convergence was fine, particularly when the exact solution was a polynomial. As shown elsewhere [66,77], all conventional finite elements of the Lagrange family can be produced by the Gordon-Coons interpolation, whereas all conventional finite elements of the Serendipity family are produced by the "boundary-only" Coons interpolation.
Moreover, it was found that the use of Lagrange polynomials, again directly referring to the nodal values u, is equivalent with the use of Bézier (Bernstein) polynomials [111]. Of course, since modern B-splines based on DeBoor's formula are by-definition equivalent with those of Schoenberg [28], all work done so far by the NTUA's group could be repeated in terms of univariateDeBoor's B-splines, in which the matrix inversion is bypassed and also local support is achieved. We stress reader's attention to the fact that in this paragraph we referred to univariate B-splines interpolation along the four boundaries of the quadrilateral patch ABCD, as well as any interboundaries if present.
The advantage of using Gordon-Coons macroelements with internal nodes with respect to conventional Lagrangian type elements, is related to the capability of dealing with patches having a different discretization to opposite sides of the patch. An illustrative example is shown in Fig. 8d (note that the opposite sides AB and CD consist of 8 and 16 subdivisions, respectively, whereas the internal nodes form their mean average, equal to 12). Therefore, it is not necessary to uniformly divide each side whereas the internal nodes are not obliged all to possess the same position in horizontal and vertical lines. For other examples we refer to [37,72,80]. A second advantage is the capability of Gordon-Coons interpolation to deal with dissimilar patches, which can be easily merged [37, pp. 330-331; 88]. Obviously, three-dimensional problems are analogous with two-dimensional ones.
In this review paper we also extended the Gordon-Coons models previously presented in [80]. We have revealed that "playing" with both the basis and the blending functions we can develop a large number of alternative Coons-elements, some of which (not all of them) have been tested in [80]. It is quite interesting that even the classical four-node (bilinear) element can be produced using hat functions for both the basis and the blending functions.
Independently of the model choice, global interpolation causes a high CPU-time, sometimes higher than the usual finite element solution [68,150]. The latter has kept researchers busy for a long time. Provatidis [73, p. 6704] was the first who realized the need for replacing the domain integration of the Galerkin-Ritz by a Global Collocation scheme where the partial differential equation is satisfied at Gauss points, taken for the entire macroelement. In this way the overall CPU-time is highly reduced.
While Galerkin-Ritz can work to complex regions, such as L-shaped or even Pi-shaped, in conjunction with piecewise linear interpolation, the latter is not possible for the GCM (it works but leads to very poor results). Lagrange polynomials and B-splines interpolation do not work along an L-shaped line because they lead to an interpolant that extremely deviates from the actual line. Therefore, the domain decomposition becomes quite necessary.
Remaining by the Global Collocation Method (GCM), our numerical experience is that the domain decomposition and the necessary coupling between adjacent macroelements do not immediately lead to acceptable results, unless an adequate number of collocation per direction is considered. A general characteristic of the Global Collocation Method is that the degrees of freedom possessing a Neumann boundary condition, or is at the interface between two adjacent macroelements, are finally eliminated and thus decrease the dimension of the system that has to be solved.
There is no doubt that the most efficient way to accurately approximate the geometry variations of an object or a mechanical component is NURBS. For example, in order to approximate the surface of an ideal sphere, it becomes necessary to divide it into at least six equal Coons patches. Similar difficulties have forced some companies to abandon Coons interpolation and replace it by NURBS [161].

Conclusions
It was shown that all the well-known interpolations that are used for the analytical description of CAD surfaces or volumes, can become the "vehicle" to derive global shape or basis functions and thus to perform engineering analysis. The only restriction is that the subregions in which the entire structure has to be subdivided should be of regular shape that can be mapped to a square, a triangle, a cube, and so on. Particularly in dynamic analysis, since the approximation of the variable is continuous in a large area or volume, the produced modes and the relevant natural frequencies are of excellent quality. The advantage of the CAD-based approach is that the degrees of freedom are associated to points that are quite necessary to determine the geometry, a fact that enables shape optimization to be efficiently performed. Concluding, all CAD-based global approximations can produce the basis functions that can be involved in any usual computational method, such as the finite element method, the boundary element method and several global collocation schemes.

Acknowledgements
This work was partially supported by the Committee of Basic Research at NTUA, Greece, through the Project: Protagoras: ΔΙΔΥΜΟ (ΔΙΕΡΕΥΝΗΣΗ ΤΗΣ ΔΥΝΑΤΟΤΗΤΑΣ ΣΥΖΕΥΞΗΣ ΤΩΝ ΜΕΘΟΔΟΛΟΓΙΩΝ CAD/CAE), EDEIL #65/1388 (Investigation on CAD/CAE integration, 1 June 2004 -28 February 2007). The proposal was submitted in May 2003. In addition, I acknowledge the PhD students who tried to deal with this matter since early 1980s at NTUA, as well as the undergraduate students who conducting their Diploma Work or MSc Thesis on this topic under my supervision. I have to point out that, primarily due to insufficient funding, and secondarily because this subject requires full concentration and dedication as well as excellent programming skills, some doctoral students discontinued and therefore the progress was slow.

Appendix A
Barnhill's interpolation to derive the linear shape functions of triangular element Let us consider the triangular element ABC≡123, as shown in Fig. 8a. Assuming the variation of the variable, F, along every side to be linear, the boundary terms shown in Fig. 4 can be written as:

Appendix B
Barnhill's interpolation to derive the quadratic shape functions of triangular element Let us consider the triangular element ABC≡123456, as shown in Fig. 8b. Since the variation of the variable, F, along every side, S (= AB,BC, CA), is quadratic, we introduce the Lagrange polynomials:    It is evident that ½ of the functions that multiply the nodal values 6 1 2 3 4 5 , , , , ,F F F F F F in (B-6) are the shape functions of the six-node finite element.
Substituting (B-1) in (B-6), the shape function ( ) 1 , , N u v w becomes: it is immediately concluded that (B-7) until (B-12) coincide with those mentioned in finite element textbooks, for example [5,6 ]. Remark: Barnhill's formula allows for the construction of arbitrary-noded triangular macroelements. In such as case, Equation (B-6) must be properly extended and then is easily programmed.