- Email: [email protected]
Spectral approximations to PDE surfaces
Computer-AidedD&gn,
Vol. 28, No. 2, pp. 145-152, 1996 Copyright 0 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 00104465/96 $15.00 + 0.00
ELSEVTER
Spectral approximations surfaces
to PDE
Malcolm I G Bloor and Michael J Wilson tion. Previously a variety of methods have been used to generate PDE surfaces. In some circumstances, for certain boundary conditions, it is possible to obtain surface solutions in closed-form6, and computationally this makes their calculation and recalculation in response to user manipulation, very efficient. Thus, the interactive manipulation in real-time of objects constructed from such solutions is a practical proposition given a work station of moderate computational power and adequate graphics performance. It has been pointed out, however, that in the case of general boundary conditions a numerical solution to the PDE must be sought ‘s7, and a number of different techniques are available, e.g. finite element*,‘, finite difference”, collocation”~12 and spectral methods’“-15. Although some of these are very efficient as far as numerical methods go, given the present state of CPU performance of the work stations commonly available, it does take longer than a fraction of a second to generate a solution, and as far as interactive design is concerned this is not ideal. Furthermore, having to obtain a PDE solution using such numerical techniques does not allow the powerful capability of the method to be fully realized. As a means of allowing interactive design in the case of general boundary conditions, this paper describes a method whereby approximate solutions to generating PDEs which satisfy the associated boundary conditions may be efficiently calculated in an explicit form. The solutions are expressed in terms of a finite sum of analytic functions (i.e. closed-form solutions) which individually satisfy the partial differential equation but not the boundary conditions, to which is added a ‘corrector’ term. The approximate solutions discussed in this paper can be made to approach the true solution to any degree of accuracy and, besides, are valid surfaces in their own right which share many of the desirable properties of PDE surfaces. An important feature of the method, and one that distinguishes it from existing solution methods, is the presence of a corrector or ‘remainder’ term as it is referred to in this paper. The purpose of this term is to ensure that at the boundary the solution exactly matches the imposed boundary conditions. Although it is natural to describe the approximation
The PDE method generates surfaces from the solutions to elliptic partial differential equations (PDEs), where boundary conditions are used to control surface shape. This paper describes a method whereby PDE surfaces may be obtained in closed form, even for the case of general boundary condi-
tions. Furthermore, the method is fast, making possible the interactive manipulations of PDE surfaces in real-time. Keywords: fourier transforms, faces, real-time
PDE surfaces, parametric sur-
INTRODUCTION The PDE method for surface generation produces surfaces as the solutions of boundary value problems, in particular as solutions to elliptic partial differential equations. The method was introduced in the context of blend generationi,’ where the object is basically to find a function which is smooth and which satisfies certain continuity conditions along its boundaries. The use of the method to generate sculptured surfaces was then discussed, and it was demonstrated that surfaces serving a wide variety of functions could be produced, e.g. ship hulls, propeller blades3,4. The PDE method is primarily a technique for surface generation rather than surface representation, and so the ultimate criteria by which its results are judged should be based upon how well the objects so described meet their objectives, whether they be aesthetic or functiona13. This paper is concerned with how PDE surfaces may be represented mathematically in order to exactly define their shape and other surface properties at any point. As mentioned above, boundary conditions imposed around the edges of a surface patch control the internal shape of the surface. The question of how the boundary conditions should be chosen has been addressed in earlier work, e.g. Reference 5, and the prime concern of this paper is the fast solution in closed form of the generating PDEs, which is important not only if the method is used to design interactively in real-time, but also as it removes any requirement for interpolaDepartment of Applied Mathematical Studies, The University of Leeds, Leeds LS2 9JT, UK Paper receiwd: 16 March 1995. Retied:
I6 June 1995
145
Spectral approximation
to PDE surfaces: M I G Bloor and M J Wilson
as a separable solution in truncated series form and a corrector term to handle the boundary conditions, the solution can be regarded as a spectral approximation’3-‘5 to the solution of the chosen PDE, in that it represents the solution in terms of globally defined functions that are infinitely differentiable. The method uses a somewhat unconventional basis for the representation of the surfaces it produces, in that it makes use of an orthogonal system of analytic functions, i.e. the Fourier system, although the use of orthogonal systems in the CADcontext is not unknown’“. Thus, for the purposes of data transfer to other CAD systems, we envisage the method being used in conjunction with a (slower) method which can calculate PDE surfaces in terms of polynomial spline functionss,i2; the process of shape design being carried out using the present method, then the surface being calculated as a collection of B-splines for transfer to another system (if necessary) once the design parameters have been chosen.
THE PDE METHOD The PDE method produces a surface as a solution to a suitably chosen partial differential equation. Past work has concentrated upon solutions to the following equation, based upon the biharmonic equation
d’
(
d2
2
X=0
Gfa2j-g
(1)
1
This equation is solved over some finite region 1IZof the (u, v) parameter plane, subject to boundary conditions on the solution X(u, U) which usually specify how X and its normal derivative dX/drn vary along c?R. The three components of the function X are the Euclidean coordinate functions of points on the surface given parameterically in terms of the two surface parameters u and u which define a coordinate system on the surface. Note that in the simplest case Equation 1 is solved independently for the x, y and z coordinates. The partial differential operator in Equation 1 represents a smoothing process in which the value of the function at any point on the surface is, in a certain sense, a weighted average of the surrounding values17. In this way a surface is obtained as a smooth transition between the boundary conditions. The parameter a controls the relative rates of smoothing between the u and u parameter directions, and for this reason has been called the smoothing parameter.
AN EFFICIENT SOLUTION METHOD We will consider the solution of Equation 1 over the (u, U) region (n: [O,11x [0,27rl, where u varies 0 -+ 27r, subject to periodic boundary conditions in the u direction, i.e. topologically, the surface is like a closed ‘band’ with the u = 0 and u = 1 isolines forming the boundary curves for the patch. If we assume that the boundary conditions take the form X(0, U) = f,(u),
X(1, U>= f,(u)
X,(0, u) = so(u), X,(1, L:) = s&L’)
146
(2) (3)
where the functions fat u), f,( c>, s&u:) and s,(u) are specified, then, by using the method of separation of variables17, the solution of Equation 1 may be written X(U,L’) =A&)
R; + c [A,Cu)cos(nu) n=l
+B,(u)sin(nu>l
(41
where the ‘coefficient’ functions A,(u) and B,(u) are of the form A,,(u) = a,,, + a,,,u + a02u2 + a03u3 A,,(u) = an,eanU + an2ueunu + a,3em”“” + a,4ue
(5) *nu
(6) B,(u) = b,,,eUnn + b,*ueanfr f bn3e-un’c + b,4ue
“I’
(7) and a,,, an2, an3, an4, h,,, bn2, hn3, h,, are constant vectors. Note that each of the three types of term, A,(u), A,(u)cos(nu) and B,(u)sin(nu) in Equation 4 satisfy Equation 1. For a given set of boundary conditions, in order to determine the various constants in the solution, it is necessary (in general) to Fourier analyse the boundary conditions and identify the various Fourier coefficients with the values of the coefficient functions A,(u) and B,(u) (and their derivatives with respect to u) at u = 0 and u = 1. In the case that the functions f,(u), f,(u), s,J 1’) and s,(u) are given (or at least expressible) as finite Fourier series, the series solution in Equation 4 is also finite, and, given the speed of fast-Fourier transforms”, this is an efficient way of obtaining the solution computationally, However, in the case that the Boundary conditions 2 and 3 are not expressible as finite Fourier series, then a solution of the form of Equation 4, now an infinite (N = a> series, is of little practical use, and other solution techniques should be used as outlined in the introduction. However, consider the form of the solution given by Equation 4. Note the coefficient functions of each Fourier term (except the zerothl may be written in the form Cd + eu)e aMu- 1) + (f+ gu)e-anw-l,
(8)
It can be seen that the first and second terms decay exponentially with u away from the boundary u = 1, while the third and fourth terms decay exponentially with u away from the boundary u = 0. Furthermore, the higher the Fourier mode, in particular the larger the value of an, the more rapid the decay of these terms as one move into the interior of the patch. Thus, Fourier modes above a certain frequency make a negligible contribution to the bulk of the surface. They are important, however, close to the boundaries where the surface patch matches up to the boundary conditions. At just what distance from a boundary the high frequency becomes negligible, depends on the value of a and the amplitudes of the high frequencies relative to the low frequencies in the boundary conditions. The observations of the above paragraph are the basis of the solution technique that is the subject of
Spectral approximation to PDE surfaces: M I G Bloor and M J Wilson
this paper. The basic idea is to approximate the solution to Equation 1 by an expression of the form
ds,( u) = R,(l, u)
N
Xh,u)
= A,(u) + c
[A,(u) cos(nu)
n=l
+B,(u)sin(nu)]
+ R(u, u)
(9)
where N is finite and in the typical problem is not large at all, e.g. considering 5 Fourier modes is often more than adequate; the function R(u, u) is a ‘remainder’ term, which is described below. The coefficient functions A,(u) and B,(u) are given by Equations 6 and 7, and may be determined from the amplitudes of the first N Fourier modes in the boundary conditions. For example, a Fourier analysis of the boundary conditions allows f, to be expressed in the form f,(u) = a,, + 5 [a, cos(nu> + b, sin(
(10)
n=l
Hence, by comparison of Equations 9 and 10, we can deduce A,(O) = a,, B,(O) = b,, etc. Thus, there is a simple identification between the values of the coefficient functions of Equation 9 and the coefficients of the Fourier series expansion of the Boundary conditions 2 and 3. The remainder function R(u, u> represents the contribution of high frequency modes to the surface: a contribution which is negligible over most of the patch if N is large enough. In this paper R(u, v> is chosen to be of the form R(u, u) = rl( u)e”” + r2( u)ue”” + r3( u)e-wu + r4( u)ue-m*
(11)
where the coefficient functions rl(u), r2(u), r3(u), r&u> are determined in the following manner. For the sake of brevity define a function F(u, u) thus F(u,u) = A,(u) + t
[A,(u)cos(nu>
n=l +
B,(u) sin(nu>]
(12)
where the A,(u) and B,(u) are obtained from a Fourier analysis of the boundary conditions as outlined above. Then define four functions df,(u), df,(u), ds,(u), ds,(u), giving the difference between the original boundary conditions and the boundary conditions satisfied by F(u, u>, thus df&)
= f,(u) - F(0, u>, df,(u) = f,(u) - F(1, u> (13)
ds,,( u> = s,-,(u> - F,(O, u), ds,( u> = sl( u> - I?,& u> (14) where the subscript u denotes partial differentiation with respect to U. Then the four functions rl(u), r&u), r&u), r4(u) are determined from the following four relations: df,( u) = R(0, u> df,(rA = R(1, u)
ds,(u) = R,(f), u> (15)
Note that the function R(u, u) is effectively expressed in terms of the boundary conditions, and by making this choice of R(u, u) we obtain an approximate solution that exactly satisfies the original boundary conditions. In choosing R(u, u) to be the form of Equation 11 we are, in effect, solving Equation 1 with a = a(n). To see this, consider the form of Solution 4 and note that each coefficient function A,(u) satisfies the ordinary differential equation
($
-.‘n2rA&,
=0
The product an in the above differential operator is the coefficient of u in the exponential terms of A,(u). We can write the solution to Equation 1 in the following form, X(u, u> = A,(u) + f [A,(u) co&u) il=l +B,(u)
sin(nu>] +
i n=N+
+
B,(u) sin(nu>]
[A,(u) cos(nu> 1
(16)
and if this is consistent with Equations 9 and 11, we are assuming that a is given by a (constant) w/n
n IN, n>N.
To ensure that the approximate Solution (9) is close to the exact solution corresponding to the case a = constant for all n, we must choose appropriately the constant w that controls the rate at which R(u, u) decays away from the boundaries. One choice is to put w = u(N + 1) which, if the (N + 1) mode is the dominant among those ignored, is close to the actual decay rate of the difference [X(u, u) - F’(u, u)]. However, there is nothing preferred about this choice for w, other values are possible and the capacity to vary w means that one can have a smoothing rate for longlength scale features in the surface that is different from the smoothing rate for small-scale features. Indeed, there is nothing to prevent one varying not only w, but also the smoothing rate for all the modes that appear explicitly in F(u, u). The choice of N will obviously affect how good an approximation Solution 9 is to the solution of the original partial differential Equation 1, and whether or not a particular value of N represents a good choice in this sense depends not only upon the boundary conditions, but ultimately on the suitability of the generated surface. However, noting that the amplitude of the nth mode decays by a factor e-l over a length scale O[(l/un)], and assuming that high frequency modes are not strongly represented in the boundary condition, in practice it seems that it is not unreasonable to 147
Spectral approximation
to PDE surfaces: M
I G Boor and M J Wilson
expect that for values of u = 1. N need not be much greater than 5. The surfaces generated by this method share many of the useful properties of the exact solution, including the fact that they both satisfy the same set of boundary conditions. Furthermore, from the point of view of visual appearance they are virtually indistinguisable and may be calculated very rapidly, which is important when designing in real-time. Having once settled on a set of boundary conditions that produces a satisfactory surface, one can then recalculate more accurately the surface using this or another method, e.g. in terms of B-splines if data transfer is an issuex.
the choice of the derivative boundary conditions; the procedure adopted here is basically that described in Reference 5 in which the u-isolines of the PDE surface are orthogonal to the boundary trimline. An example of such a blend is shown in Figure I. Note that these results were prepared using software which allows the blend surface to be changed interactively in real-time in response to changes in the boundary conditions and derivatives. Due to the finite size of the facets into which the surfaces have been broken for the purposes of rendering, the edges of the blend look slightly ‘ragged’, but in reality are perfectly smooth.
RESULTS
Design of a generic aircraft shape
We will now consider two examples of the use of this method. The first is concerned with the generation of a secondary blending surface between two primary Bspline surfaces, and the second is concerned with the use of the method in parametrizing aircraft geometries.
In this example we consider the design of a generic aircraft shape which is made up of three patches: a fuselage, an inner wing, and an outer wing. We consider in particular the double-delta configuration characteristic of a supersonic aircraft. Although this work has been described elsewhere”, we shall give a fair amount of detail here to illustrate the way in which the boundary conditions are chosen in order to join multiple PDE surfaces to make complex geometries. For simplicity we will use a fuselage that is defined algebraically (see below). The inner and outer wings will be generated using the PDE method. The ‘characterlines’ which form the boundaries between adjacent surface patches are:
Blending between two B-spline surfaces The following describes the way in which the method described in this paper may be used to generate a blend between two B-spline surfaces. The details of this application, in particular the calculation of the function and boundary conditions, have been described elsewhere5-albeit in the context of blending between quadric surfaces-and the interested reader is directed to this reference; here we will merely give an outline of the method. We shall take the B-spline primary surfaces as given, and concentrate upon finding a blending surface between them that joins with continuity of surface normal. To generate the blend surface, we define ‘trimlines’ on the surface of each B-spline surface. This is done by means of a 2D B-spline curve defined in the parameter space of each primary surface. For instance. if XA(r3,~> is one of the B-spline surfaces, then a trimline XA[0(f),4(t>] may be defined on it in terms of a curve parameter t, where the functions 8 = 8 (t 1, I#J= 4(f) relating the surface parameters to the curve parameter are the usual sum over a set of B-spline basis functions B#), thus [e(f),c#df)l=
&;B,W i
where the pi are control points in the (e,+I plane. A similar procedure is then followed for the second Bspline surface XB. Now, since we have an exact representation of each primary surface and of each trimline in terms of Bsplines, we have the function boundary conditions defined at every point on the trimlines. We can also obtain derivative conditions at each point on the trimlines by ensuring that the surface normal to the PDE blend (determined by X, X X,) is parallel to the surface normal of the B-spline surface along the trimline. Note that within the constraints set by continuity of surface normal, there is still considerable latitude in 148
( I ) the curve where the inner and outer wing meet, (2) the curve where the inner wing meets the fuselage, and (3) a curve at the tip of the outer wing.
Function boundary conditions The fuselage is generated as a surface of revolution whose axis is parallel to the x axis, and where the y and z coordinates of points on the surface are related by \”
+
z? = a:
(17)
with
u( x) =aosin
18 Z(17*+ L
1) +a,sin I
37* i
+ 1)
1
(18)
where a,,, a, are constants and x is a parameter which lies in the range 0 5 x I 1. Note that as x varies in the range 1 + 0 we move from the front towards the rear of the aircraft, and the cylindrically symmetric fuselage exhibits the waisted profile characteristic of an aircraft designed for supersonic flight. The curve where the outer wing and inner wing meet we will take to be a plane curve (z = constant) having the shape of a simple airfoil, described parametrically
Spectral approximation
to PDE surfaces: M
I G Boor and M J Wilson
Fi@Ire 1 PDE blend between B-spline primary surfaces
thu
+ (6.75)(cam)x(8)
X
>=ch*[l-cos(8)1/2
Y
t ) = - -sin(e)
2
x[ch
-x(O)l[ch -xuN/ch3
z=a,+h,
+ E * sin(28) where the parameter
(19) 8 varies in the range 0 I 0 I : 2lr, 149
Spectral approximation to PDE surfaces: M I G Bloor and M J Wilson cam, ch, hi, and E are parameters controlling camber, chord length, inner section span and the rounding of the trailing edge. The second characterline lies on the surface of the fuselage, and its projection onto a vertical plane containing the fuselage’s axis is an airfoil of the same type as characterline 1, but scaled by a factor (b/c/r) and offset with respect to first characterline by the vector (Gxd, yd). It is given parametrically by the equations
yf(e)
=yd+
z,(O) =
&Y(B)
I/u’r x(e)1
-Yj(@
(20)
where a is given in terms of x by Equation 19, and x is given in terms of 8 by bx(O)
x(e) = -
rl ch
boundary conditions on X,L(= d X/a u). The form of these derivative boundary conditions is chosen so that there is continuity of surface normal between adjacent surface patches. Now, the surface normal is determined by the vector product X, X X,. and, in the example, X, is given on the boundary by the conditions on X, which means that the surface normal at the edge of a surface patch is determined by the boundary conditions we choose to specify for X,,. An example of the aircraft shapes that may readily be produced is shown in Figure 2. It has been produced by a program that allows a designer, sitting at a work station, to interactively define a wide range of shapes in real-time. Throughout the geometry changes, the continuity of the surfaces patches that define the aircraft shape is automatically maintained. The parametrization of the shape gives rise to a surface mesh that may readily be used as the basis for physical analysis and also for the manufacture of models for physical testing, using such rapid proto techniques as stereolithography or laser sintering %P’ .
xte
(21)
+rl
CONCLUSION where b, rl, xte and rl are constants. The third characterline lies at the tip of the outer wing. It is given parametrically by the equations
rxtlx(6,)
1
x,(e) =xt + ~chJCOshd xtlx(
y,(B)=-
---$y
[ z,=a,+h,
+h,
8> 1
sin(ao) (22)
We can see from the above equations that as 0 varies in the range 0 * 27r, the wingtip is a closed straight line of length xtl. For simplicity we will assume that on the boundaries of the surface patches the curve parameter and the surface parameter are identical, i.e. 8 = L’,so that for the inner wing X(0, ~1) is given by Equation 19, and X(1, u) is given by Equations 20; similarly for the outer wing. We will consider wings that are closed at the trailing edge and hence we look for solutions with the property that X(u, 0) = X(u, 2~). Two wing surfaces are generated between these three characterlines: the inner and outer surfaces of the double-delta. The inner wing is generated by solving Equation 1 using boundary conditions obtained from the Characterlines 19 and 20, and the outer wing is generated by solving Equation 1 using boundary conditions obtained from the Characterlines 20 and 22. The latter can be represented by a finite Fourier series, whereas the former requires the approach outlined in this paper essentially due to Boundary Conditions 20. Derivative boundary conditions solution
for PDE
Equation 1 is 4th order, hence we require boundary conditions on the normal derivatives of X(u, c) in the (u, u) parameter plane, which is this case means 150
We have described a method whereby approximate solutions to PDE surfaces can be calculated. It involves a Fourier analysis of the boundary conditions which, using fast Fourier transform techniques”, can be done very efficiently. The solution is expressed in terms of a finite number of analytic functions-which are closedform solutions to the partial differential equation though not the boundary conditions-plus a remainder function (which itself will be analytic if the boundary conditions are expressed in terms of analytic functions). The remainder function is included in the solution to ensure that, although approximate in the surface interior, it satisfies the boundary conditions exactly (in whatever form they are expressed). Thus we have an expression for the surface which gives all points on it as explicit functions of the surface parameters u and ~1. This facilitates the physical analysis of the object and its manufacture by both conventional techniques and by novel layering techniques2”. Using this technique, we can obtain solutions for PDE surfaces with general boundary conditions rapidl? enough to be able to design interactively in real-time . For example, the method described in this paper has been implemented in software that can run on Silicon Graphics work stations, and the results given in this paper have been created using this software. The PDE surfaces change at the same rate as any alterations in the design parameters effected via the user interface, i.e. at moving-picture speed, since the solution for 1000 facets is of the order of 0.01s on a Silicon Graphics 4D420. This is in constrast to other methods that have been used to calculate PDE surfaces. For instance, if data transfer to other CAD systems is an issue, it is possible to represent PDE surfaces in terms of B-splines using the finite-element method’. Using this approach, calculating a PDE surface solution may take, typically, from the order of a second to a minute, depending on the size of the B-spline control mesh. The method can be viewed as belonging to the class of spectral methods for the solution of partial differen-
Spectral approximation to PDE surfaces: M
Figure 2
PDE aircraft
I G Boor and M J Wilson
geometry
tial equations’3%1s where the expansion functions are based upon the Fourier system. Such methods are often used for the modelling of physical systems, e.g. fluid flow14 but unlike these applications the geometric context in which the present method will be used is
unlikely to involve the presence of discontinuities in the function boundary conditions, i.e. the boundaries of the PDE surfaces will be continuous curves. Thus, the problem of the Gibbs phenomenon is unlikely to arise from discontinuities in the function boundary condi151
Spectral approximation
to PDE surfaces: M I G Bloor and M J Wilson
tions. This paper has only considered periodic solutions of the PDE where the conditions have contained no singularities. However, there is the possibility of discontinuities arising (fairly) naturally in the derivative boundary conditions when considering boundary curves with corners. In a later paper it will be shown how the method can be adapted to produce PDE surfaces having boundary curves with corners (or other types of singularity), so that a wide range of complicated geometries can be produced, e.g. inlet ports to internal combustion engines*i .
12 13
14 15
lb 17 1x
ACKNOWLEDGEMENTS The authors would like to thank Dr RE Smith of the NASA Langley Research Center for his encouragement in the area of aircraft design which was supported by NASA grant NAGW-3198.
19 20
21
Bloor, M 1 G and Wilson, M J ‘Representing PDE surfaces m terms of B-splines’ Comput.-Aided Des. Vol22 (19901 pp 324-331 Gottleib, D and Ortega, S ‘Numerical analysis of spectral methods: theory and applications’ CBMS-NSF Regional Conf Series in Applied Mathematics No. 26 SIAM, Philadelphia (19771 Canuto, C, Hussaini, M Y, Quateroni, A and Zang, T Spectrul Methods in Fluid Dynamics Springer, Berlin (1988) Quateroni, A and Valli, A Numerical Approximations of furtial Differential Equations Springer Series in Computational Mathematics 23, Berlin (19941 Lachance, M A ‘Chebyschev errors for parametric surfaces’ Comput. Aided Geom. Des. Vol 5 (1988) pp 195-208 Zauderer, E Partial Differential Equations of Applied Mathematics Wiley Interscience, New York (1983) Elliot, D F and Rao, K R Fast Transforms Academic Press. New York (1982) Bloor, M I G and Wilson, M J ‘The efficient parametrization of generic aircraft geometries’ J Aircraft, to appear (1995) Kai, C C ‘Three-dimensional rapid prototyping technologies and key development areas’ Comput. Control Engng .I. Vol S No 4 (1994) pp 200-206 Bloor, M I G and Wilson, M J ‘Complex PDE surface gcneration for analysis and manufacture’ Computing Vol 10 (199% pp 61-77
REFERENCES I
2 3
4
5
6
7
8
9 IO 11
152
Bloor, M I G and Wilson, M J ‘Generating blend surfaces using partial differential equations’ Comput. Aided. Des. Vol 21 No 3 (19891 pp 165-171 Hoschek, J and Lasser, D Fudamentais of Computer Aided Geometric Design A K Peters, Wellesley, MA (1993) Lowe, T W, Bloor, M I G and Wilson, M J, ‘The automatic design of hull surface geometries’ J. Ship Res. Vol 38 No 4 (1994) pp 319-328 Dekanski, C, Bloor, M I G and Wilson, M J ‘The representation of marine propeller blades using the PDE method’ .I. Ship Res. to appear (1995) Bloor, M I G and Wilson, M J, ‘Interactive design using partial differential equations’ in Sapidis, N (Ed.1 Designing Fair Curces and Surfaces SIAM, Philadelphia (19941 pp 231-25 1 Bloor, M I G and Wilson, M J ‘Using partial differential equations to generate free-form surfaces’ Cornput-Aided Des. Vol 22 (1990) pp 202-212 Vida, J, Martin, R R and Varady, T ‘A survey of blending methods that use parametric surfaces’ Cornput-Aided Des. Vol 26 No 5 (19941 pp 341-365 Brown, J M, Bloor, M I G, Bloor, M S and Wilson, M J, ‘Generation and modification of non-uniform B-spline surface approximations to PDE surfaces using the finite-element method, in Ravani, B (Ed.) Advances in Design Automation Vol 1 Computer Aided and Computational Design ASME (1990) pp 265-272 Prenter, Splines and VariationalMethods Wiley-Interscience, New York (1975) Burden, R L and Faires, J D Numerical Analysis (3rd edn) Prindle, Weber and Schmidt, Boston (19811 Ortega, J M and Poole, W G An Introduction to Numerical Methods for Differential Equations Academic Press, New York (1981)
Malcolm Bloor graduated with a BSc in mathematics from Manchester Unitersity in 1962. He obtained a PhD in applied mathematics from Manchester Uniwrsity in 1966. He was appointed as an assistant lecturer in the Department of Applied Mathematical Studies in 1964. He cuwently has a Chair in Mathematical Engineering at Leeds His research interests include computeraided design, fluid mechanics, industrial 1 mathematics, and biological mathemattcs
Michael Wilson graduated in 1980 with a BA in natural science from Cambridge University After research as pun of the Radio Astronomy Group in the Caoendish Laboratory, he obtained his PhD in I984 from Cambridge Unitersity. He joined the Department of Applied Mathematical Studies at Leeds University as a research fellow in 1983, and was appointed to a lecturing position in the same department in 1986. He is currently a senior lecturer. His research interests include computer-aided design, computational fluid dynamics, industrial mathematics and astrophysics.
Vol. 28, No. 2, pp. 145-152, 1996 Copyright 0 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 00104465/96 $15.00 + 0.00
ELSEVTER
Spectral approximations surfaces
to PDE
Malcolm I G Bloor and Michael J Wilson tion. Previously a variety of methods have been used to generate PDE surfaces. In some circumstances, for certain boundary conditions, it is possible to obtain surface solutions in closed-form6, and computationally this makes their calculation and recalculation in response to user manipulation, very efficient. Thus, the interactive manipulation in real-time of objects constructed from such solutions is a practical proposition given a work station of moderate computational power and adequate graphics performance. It has been pointed out, however, that in the case of general boundary conditions a numerical solution to the PDE must be sought ‘s7, and a number of different techniques are available, e.g. finite element*,‘, finite difference”, collocation”~12 and spectral methods’“-15. Although some of these are very efficient as far as numerical methods go, given the present state of CPU performance of the work stations commonly available, it does take longer than a fraction of a second to generate a solution, and as far as interactive design is concerned this is not ideal. Furthermore, having to obtain a PDE solution using such numerical techniques does not allow the powerful capability of the method to be fully realized. As a means of allowing interactive design in the case of general boundary conditions, this paper describes a method whereby approximate solutions to generating PDEs which satisfy the associated boundary conditions may be efficiently calculated in an explicit form. The solutions are expressed in terms of a finite sum of analytic functions (i.e. closed-form solutions) which individually satisfy the partial differential equation but not the boundary conditions, to which is added a ‘corrector’ term. The approximate solutions discussed in this paper can be made to approach the true solution to any degree of accuracy and, besides, are valid surfaces in their own right which share many of the desirable properties of PDE surfaces. An important feature of the method, and one that distinguishes it from existing solution methods, is the presence of a corrector or ‘remainder’ term as it is referred to in this paper. The purpose of this term is to ensure that at the boundary the solution exactly matches the imposed boundary conditions. Although it is natural to describe the approximation
The PDE method generates surfaces from the solutions to elliptic partial differential equations (PDEs), where boundary conditions are used to control surface shape. This paper describes a method whereby PDE surfaces may be obtained in closed form, even for the case of general boundary condi-
tions. Furthermore, the method is fast, making possible the interactive manipulations of PDE surfaces in real-time. Keywords: fourier transforms, faces, real-time
PDE surfaces, parametric sur-
INTRODUCTION The PDE method for surface generation produces surfaces as the solutions of boundary value problems, in particular as solutions to elliptic partial differential equations. The method was introduced in the context of blend generationi,’ where the object is basically to find a function which is smooth and which satisfies certain continuity conditions along its boundaries. The use of the method to generate sculptured surfaces was then discussed, and it was demonstrated that surfaces serving a wide variety of functions could be produced, e.g. ship hulls, propeller blades3,4. The PDE method is primarily a technique for surface generation rather than surface representation, and so the ultimate criteria by which its results are judged should be based upon how well the objects so described meet their objectives, whether they be aesthetic or functiona13. This paper is concerned with how PDE surfaces may be represented mathematically in order to exactly define their shape and other surface properties at any point. As mentioned above, boundary conditions imposed around the edges of a surface patch control the internal shape of the surface. The question of how the boundary conditions should be chosen has been addressed in earlier work, e.g. Reference 5, and the prime concern of this paper is the fast solution in closed form of the generating PDEs, which is important not only if the method is used to design interactively in real-time, but also as it removes any requirement for interpolaDepartment of Applied Mathematical Studies, The University of Leeds, Leeds LS2 9JT, UK Paper receiwd: 16 March 1995. Retied:
I6 June 1995
145
Spectral approximation
to PDE surfaces: M I G Bloor and M J Wilson
as a separable solution in truncated series form and a corrector term to handle the boundary conditions, the solution can be regarded as a spectral approximation’3-‘5 to the solution of the chosen PDE, in that it represents the solution in terms of globally defined functions that are infinitely differentiable. The method uses a somewhat unconventional basis for the representation of the surfaces it produces, in that it makes use of an orthogonal system of analytic functions, i.e. the Fourier system, although the use of orthogonal systems in the CADcontext is not unknown’“. Thus, for the purposes of data transfer to other CAD systems, we envisage the method being used in conjunction with a (slower) method which can calculate PDE surfaces in terms of polynomial spline functionss,i2; the process of shape design being carried out using the present method, then the surface being calculated as a collection of B-splines for transfer to another system (if necessary) once the design parameters have been chosen.
THE PDE METHOD The PDE method produces a surface as a solution to a suitably chosen partial differential equation. Past work has concentrated upon solutions to the following equation, based upon the biharmonic equation
d’
(
d2
2
X=0
Gfa2j-g
(1)
1
This equation is solved over some finite region 1IZof the (u, v) parameter plane, subject to boundary conditions on the solution X(u, U) which usually specify how X and its normal derivative dX/drn vary along c?R. The three components of the function X are the Euclidean coordinate functions of points on the surface given parameterically in terms of the two surface parameters u and u which define a coordinate system on the surface. Note that in the simplest case Equation 1 is solved independently for the x, y and z coordinates. The partial differential operator in Equation 1 represents a smoothing process in which the value of the function at any point on the surface is, in a certain sense, a weighted average of the surrounding values17. In this way a surface is obtained as a smooth transition between the boundary conditions. The parameter a controls the relative rates of smoothing between the u and u parameter directions, and for this reason has been called the smoothing parameter.
AN EFFICIENT SOLUTION METHOD We will consider the solution of Equation 1 over the (u, U) region (n: [O,11x [0,27rl, where u varies 0 -+ 27r, subject to periodic boundary conditions in the u direction, i.e. topologically, the surface is like a closed ‘band’ with the u = 0 and u = 1 isolines forming the boundary curves for the patch. If we assume that the boundary conditions take the form X(0, U) = f,(u),
X(1, U>= f,(u)
X,(0, u) = so(u), X,(1, L:) = s&L’)
146
(2) (3)
where the functions fat u), f,( c>, s&u:) and s,(u) are specified, then, by using the method of separation of variables17, the solution of Equation 1 may be written X(U,L’) =A&)
R; + c [A,Cu)cos(nu) n=l
+B,(u)sin(nu>l
(41
where the ‘coefficient’ functions A,(u) and B,(u) are of the form A,,(u) = a,,, + a,,,u + a02u2 + a03u3 A,,(u) = an,eanU + an2ueunu + a,3em”“” + a,4ue
(5) *nu
(6) B,(u) = b,,,eUnn + b,*ueanfr f bn3e-un’c + b,4ue
“I’
(7) and a,,, an2, an3, an4, h,,, bn2, hn3, h,, are constant vectors. Note that each of the three types of term, A,(u), A,(u)cos(nu) and B,(u)sin(nu) in Equation 4 satisfy Equation 1. For a given set of boundary conditions, in order to determine the various constants in the solution, it is necessary (in general) to Fourier analyse the boundary conditions and identify the various Fourier coefficients with the values of the coefficient functions A,(u) and B,(u) (and their derivatives with respect to u) at u = 0 and u = 1. In the case that the functions f,(u), f,(u), s,J 1’) and s,(u) are given (or at least expressible) as finite Fourier series, the series solution in Equation 4 is also finite, and, given the speed of fast-Fourier transforms”, this is an efficient way of obtaining the solution computationally, However, in the case that the Boundary conditions 2 and 3 are not expressible as finite Fourier series, then a solution of the form of Equation 4, now an infinite (N = a> series, is of little practical use, and other solution techniques should be used as outlined in the introduction. However, consider the form of the solution given by Equation 4. Note the coefficient functions of each Fourier term (except the zerothl may be written in the form Cd + eu)e aMu- 1) + (f+ gu)e-anw-l,
(8)
It can be seen that the first and second terms decay exponentially with u away from the boundary u = 1, while the third and fourth terms decay exponentially with u away from the boundary u = 0. Furthermore, the higher the Fourier mode, in particular the larger the value of an, the more rapid the decay of these terms as one move into the interior of the patch. Thus, Fourier modes above a certain frequency make a negligible contribution to the bulk of the surface. They are important, however, close to the boundaries where the surface patch matches up to the boundary conditions. At just what distance from a boundary the high frequency becomes negligible, depends on the value of a and the amplitudes of the high frequencies relative to the low frequencies in the boundary conditions. The observations of the above paragraph are the basis of the solution technique that is the subject of
Spectral approximation to PDE surfaces: M I G Bloor and M J Wilson
this paper. The basic idea is to approximate the solution to Equation 1 by an expression of the form
ds,( u) = R,(l, u)
N
Xh,u)
= A,(u) + c
[A,(u) cos(nu)
n=l
+B,(u)sin(nu)]
+ R(u, u)
(9)
where N is finite and in the typical problem is not large at all, e.g. considering 5 Fourier modes is often more than adequate; the function R(u, u) is a ‘remainder’ term, which is described below. The coefficient functions A,(u) and B,(u) are given by Equations 6 and 7, and may be determined from the amplitudes of the first N Fourier modes in the boundary conditions. For example, a Fourier analysis of the boundary conditions allows f, to be expressed in the form f,(u) = a,, + 5 [a, cos(nu> + b, sin(
(10)
n=l
Hence, by comparison of Equations 9 and 10, we can deduce A,(O) = a,, B,(O) = b,, etc. Thus, there is a simple identification between the values of the coefficient functions of Equation 9 and the coefficients of the Fourier series expansion of the Boundary conditions 2 and 3. The remainder function R(u, u> represents the contribution of high frequency modes to the surface: a contribution which is negligible over most of the patch if N is large enough. In this paper R(u, v> is chosen to be of the form R(u, u) = rl( u)e”” + r2( u)ue”” + r3( u)e-wu + r4( u)ue-m*
(11)
where the coefficient functions rl(u), r2(u), r3(u), r&u> are determined in the following manner. For the sake of brevity define a function F(u, u) thus F(u,u) = A,(u) + t
[A,(u)cos(nu>
n=l +
B,(u) sin(nu>]
(12)
where the A,(u) and B,(u) are obtained from a Fourier analysis of the boundary conditions as outlined above. Then define four functions df,(u), df,(u), ds,(u), ds,(u), giving the difference between the original boundary conditions and the boundary conditions satisfied by F(u, u>, thus df&)
= f,(u) - F(0, u>, df,(u) = f,(u) - F(1, u> (13)
ds,,( u> = s,-,(u> - F,(O, u), ds,( u> = sl( u> - I?,& u> (14) where the subscript u denotes partial differentiation with respect to U. Then the four functions rl(u), r&u), r&u), r4(u) are determined from the following four relations: df,( u) = R(0, u> df,(rA = R(1, u)
ds,(u) = R,(f), u> (15)
Note that the function R(u, u) is effectively expressed in terms of the boundary conditions, and by making this choice of R(u, u) we obtain an approximate solution that exactly satisfies the original boundary conditions. In choosing R(u, u) to be the form of Equation 11 we are, in effect, solving Equation 1 with a = a(n). To see this, consider the form of Solution 4 and note that each coefficient function A,(u) satisfies the ordinary differential equation
($
-.‘n2rA&,
=0
The product an in the above differential operator is the coefficient of u in the exponential terms of A,(u). We can write the solution to Equation 1 in the following form, X(u, u> = A,(u) + f [A,(u) co&u) il=l +B,(u)
sin(nu>] +
i n=N+
+
B,(u) sin(nu>]
[A,(u) cos(nu> 1
(16)
and if this is consistent with Equations 9 and 11, we are assuming that a is given by a (constant) w/n
n IN, n>N.
To ensure that the approximate Solution (9) is close to the exact solution corresponding to the case a = constant for all n, we must choose appropriately the constant w that controls the rate at which R(u, u) decays away from the boundaries. One choice is to put w = u(N + 1) which, if the (N + 1) mode is the dominant among those ignored, is close to the actual decay rate of the difference [X(u, u) - F’(u, u)]. However, there is nothing preferred about this choice for w, other values are possible and the capacity to vary w means that one can have a smoothing rate for longlength scale features in the surface that is different from the smoothing rate for small-scale features. Indeed, there is nothing to prevent one varying not only w, but also the smoothing rate for all the modes that appear explicitly in F(u, u). The choice of N will obviously affect how good an approximation Solution 9 is to the solution of the original partial differential Equation 1, and whether or not a particular value of N represents a good choice in this sense depends not only upon the boundary conditions, but ultimately on the suitability of the generated surface. However, noting that the amplitude of the nth mode decays by a factor e-l over a length scale O[(l/un)], and assuming that high frequency modes are not strongly represented in the boundary condition, in practice it seems that it is not unreasonable to 147
Spectral approximation
to PDE surfaces: M
I G Boor and M J Wilson
expect that for values of u = 1. N need not be much greater than 5. The surfaces generated by this method share many of the useful properties of the exact solution, including the fact that they both satisfy the same set of boundary conditions. Furthermore, from the point of view of visual appearance they are virtually indistinguisable and may be calculated very rapidly, which is important when designing in real-time. Having once settled on a set of boundary conditions that produces a satisfactory surface, one can then recalculate more accurately the surface using this or another method, e.g. in terms of B-splines if data transfer is an issuex.
the choice of the derivative boundary conditions; the procedure adopted here is basically that described in Reference 5 in which the u-isolines of the PDE surface are orthogonal to the boundary trimline. An example of such a blend is shown in Figure I. Note that these results were prepared using software which allows the blend surface to be changed interactively in real-time in response to changes in the boundary conditions and derivatives. Due to the finite size of the facets into which the surfaces have been broken for the purposes of rendering, the edges of the blend look slightly ‘ragged’, but in reality are perfectly smooth.
RESULTS
Design of a generic aircraft shape
We will now consider two examples of the use of this method. The first is concerned with the generation of a secondary blending surface between two primary Bspline surfaces, and the second is concerned with the use of the method in parametrizing aircraft geometries.
In this example we consider the design of a generic aircraft shape which is made up of three patches: a fuselage, an inner wing, and an outer wing. We consider in particular the double-delta configuration characteristic of a supersonic aircraft. Although this work has been described elsewhere”, we shall give a fair amount of detail here to illustrate the way in which the boundary conditions are chosen in order to join multiple PDE surfaces to make complex geometries. For simplicity we will use a fuselage that is defined algebraically (see below). The inner and outer wings will be generated using the PDE method. The ‘characterlines’ which form the boundaries between adjacent surface patches are:
Blending between two B-spline surfaces The following describes the way in which the method described in this paper may be used to generate a blend between two B-spline surfaces. The details of this application, in particular the calculation of the function and boundary conditions, have been described elsewhere5-albeit in the context of blending between quadric surfaces-and the interested reader is directed to this reference; here we will merely give an outline of the method. We shall take the B-spline primary surfaces as given, and concentrate upon finding a blending surface between them that joins with continuity of surface normal. To generate the blend surface, we define ‘trimlines’ on the surface of each B-spline surface. This is done by means of a 2D B-spline curve defined in the parameter space of each primary surface. For instance. if XA(r3,~> is one of the B-spline surfaces, then a trimline XA[0(f),4(t>] may be defined on it in terms of a curve parameter t, where the functions 8 = 8 (t 1, I#J= 4(f) relating the surface parameters to the curve parameter are the usual sum over a set of B-spline basis functions B#), thus [e(f),c#df)l=
&;B,W i
where the pi are control points in the (e,+I plane. A similar procedure is then followed for the second Bspline surface XB. Now, since we have an exact representation of each primary surface and of each trimline in terms of Bsplines, we have the function boundary conditions defined at every point on the trimlines. We can also obtain derivative conditions at each point on the trimlines by ensuring that the surface normal to the PDE blend (determined by X, X X,) is parallel to the surface normal of the B-spline surface along the trimline. Note that within the constraints set by continuity of surface normal, there is still considerable latitude in 148
( I ) the curve where the inner and outer wing meet, (2) the curve where the inner wing meets the fuselage, and (3) a curve at the tip of the outer wing.
Function boundary conditions The fuselage is generated as a surface of revolution whose axis is parallel to the x axis, and where the y and z coordinates of points on the surface are related by \”
+
z? = a:
(17)
with
u( x) =aosin
18 Z(17*+ L
1) +a,sin I
37* i
+ 1)
1
(18)
where a,,, a, are constants and x is a parameter which lies in the range 0 5 x I 1. Note that as x varies in the range 1 + 0 we move from the front towards the rear of the aircraft, and the cylindrically symmetric fuselage exhibits the waisted profile characteristic of an aircraft designed for supersonic flight. The curve where the outer wing and inner wing meet we will take to be a plane curve (z = constant) having the shape of a simple airfoil, described parametrically
Spectral approximation
to PDE surfaces: M
I G Boor and M J Wilson
Fi@Ire 1 PDE blend between B-spline primary surfaces
thu
+ (6.75)(cam)x(8)
X
>=ch*[l-cos(8)1/2
Y
t ) = - -sin(e)
2
x[ch
-x(O)l[ch -xuN/ch3
z=a,+h,
+ E * sin(28) where the parameter
(19) 8 varies in the range 0 I 0 I : 2lr, 149
Spectral approximation to PDE surfaces: M I G Bloor and M J Wilson cam, ch, hi, and E are parameters controlling camber, chord length, inner section span and the rounding of the trailing edge. The second characterline lies on the surface of the fuselage, and its projection onto a vertical plane containing the fuselage’s axis is an airfoil of the same type as characterline 1, but scaled by a factor (b/c/r) and offset with respect to first characterline by the vector (Gxd, yd). It is given parametrically by the equations
yf(e)
=yd+
z,(O) =
&Y(B)
I/u’r x(e)1
-Yj(@
(20)
where a is given in terms of x by Equation 19, and x is given in terms of 8 by bx(O)
x(e) = -
rl ch
boundary conditions on X,L(= d X/a u). The form of these derivative boundary conditions is chosen so that there is continuity of surface normal between adjacent surface patches. Now, the surface normal is determined by the vector product X, X X,. and, in the example, X, is given on the boundary by the conditions on X, which means that the surface normal at the edge of a surface patch is determined by the boundary conditions we choose to specify for X,,. An example of the aircraft shapes that may readily be produced is shown in Figure 2. It has been produced by a program that allows a designer, sitting at a work station, to interactively define a wide range of shapes in real-time. Throughout the geometry changes, the continuity of the surfaces patches that define the aircraft shape is automatically maintained. The parametrization of the shape gives rise to a surface mesh that may readily be used as the basis for physical analysis and also for the manufacture of models for physical testing, using such rapid proto techniques as stereolithography or laser sintering %P’ .
xte
(21)
+rl
CONCLUSION where b, rl, xte and rl are constants. The third characterline lies at the tip of the outer wing. It is given parametrically by the equations
rxtlx(6,)
1
x,(e) =xt + ~chJCOshd xtlx(
y,(B)=-
---$y
[ z,=a,+h,
+h,
8> 1
sin(ao) (22)
We can see from the above equations that as 0 varies in the range 0 * 27r, the wingtip is a closed straight line of length xtl. For simplicity we will assume that on the boundaries of the surface patches the curve parameter and the surface parameter are identical, i.e. 8 = L’,so that for the inner wing X(0, ~1) is given by Equation 19, and X(1, u) is given by Equations 20; similarly for the outer wing. We will consider wings that are closed at the trailing edge and hence we look for solutions with the property that X(u, 0) = X(u, 2~). Two wing surfaces are generated between these three characterlines: the inner and outer surfaces of the double-delta. The inner wing is generated by solving Equation 1 using boundary conditions obtained from the Characterlines 19 and 20, and the outer wing is generated by solving Equation 1 using boundary conditions obtained from the Characterlines 20 and 22. The latter can be represented by a finite Fourier series, whereas the former requires the approach outlined in this paper essentially due to Boundary Conditions 20. Derivative boundary conditions solution
for PDE
Equation 1 is 4th order, hence we require boundary conditions on the normal derivatives of X(u, c) in the (u, u) parameter plane, which is this case means 150
We have described a method whereby approximate solutions to PDE surfaces can be calculated. It involves a Fourier analysis of the boundary conditions which, using fast Fourier transform techniques”, can be done very efficiently. The solution is expressed in terms of a finite number of analytic functions-which are closedform solutions to the partial differential equation though not the boundary conditions-plus a remainder function (which itself will be analytic if the boundary conditions are expressed in terms of analytic functions). The remainder function is included in the solution to ensure that, although approximate in the surface interior, it satisfies the boundary conditions exactly (in whatever form they are expressed). Thus we have an expression for the surface which gives all points on it as explicit functions of the surface parameters u and ~1. This facilitates the physical analysis of the object and its manufacture by both conventional techniques and by novel layering techniques2”. Using this technique, we can obtain solutions for PDE surfaces with general boundary conditions rapidl? enough to be able to design interactively in real-time . For example, the method described in this paper has been implemented in software that can run on Silicon Graphics work stations, and the results given in this paper have been created using this software. The PDE surfaces change at the same rate as any alterations in the design parameters effected via the user interface, i.e. at moving-picture speed, since the solution for 1000 facets is of the order of 0.01s on a Silicon Graphics 4D420. This is in constrast to other methods that have been used to calculate PDE surfaces. For instance, if data transfer to other CAD systems is an issue, it is possible to represent PDE surfaces in terms of B-splines using the finite-element method’. Using this approach, calculating a PDE surface solution may take, typically, from the order of a second to a minute, depending on the size of the B-spline control mesh. The method can be viewed as belonging to the class of spectral methods for the solution of partial differen-
Spectral approximation to PDE surfaces: M
Figure 2
PDE aircraft
I G Boor and M J Wilson
geometry
tial equations’3%1s where the expansion functions are based upon the Fourier system. Such methods are often used for the modelling of physical systems, e.g. fluid flow14 but unlike these applications the geometric context in which the present method will be used is
unlikely to involve the presence of discontinuities in the function boundary conditions, i.e. the boundaries of the PDE surfaces will be continuous curves. Thus, the problem of the Gibbs phenomenon is unlikely to arise from discontinuities in the function boundary condi151
Spectral approximation
to PDE surfaces: M I G Bloor and M J Wilson
tions. This paper has only considered periodic solutions of the PDE where the conditions have contained no singularities. However, there is the possibility of discontinuities arising (fairly) naturally in the derivative boundary conditions when considering boundary curves with corners. In a later paper it will be shown how the method can be adapted to produce PDE surfaces having boundary curves with corners (or other types of singularity), so that a wide range of complicated geometries can be produced, e.g. inlet ports to internal combustion engines*i .
12 13
14 15
lb 17 1x
ACKNOWLEDGEMENTS The authors would like to thank Dr RE Smith of the NASA Langley Research Center for his encouragement in the area of aircraft design which was supported by NASA grant NAGW-3198.
19 20
21
Bloor, M 1 G and Wilson, M J ‘Representing PDE surfaces m terms of B-splines’ Comput.-Aided Des. Vol22 (19901 pp 324-331 Gottleib, D and Ortega, S ‘Numerical analysis of spectral methods: theory and applications’ CBMS-NSF Regional Conf Series in Applied Mathematics No. 26 SIAM, Philadelphia (19771 Canuto, C, Hussaini, M Y, Quateroni, A and Zang, T Spectrul Methods in Fluid Dynamics Springer, Berlin (1988) Quateroni, A and Valli, A Numerical Approximations of furtial Differential Equations Springer Series in Computational Mathematics 23, Berlin (19941 Lachance, M A ‘Chebyschev errors for parametric surfaces’ Comput. Aided Geom. Des. Vol 5 (1988) pp 195-208 Zauderer, E Partial Differential Equations of Applied Mathematics Wiley Interscience, New York (1983) Elliot, D F and Rao, K R Fast Transforms Academic Press. New York (1982) Bloor, M I G and Wilson, M J ‘The efficient parametrization of generic aircraft geometries’ J Aircraft, to appear (1995) Kai, C C ‘Three-dimensional rapid prototyping technologies and key development areas’ Comput. Control Engng .I. Vol S No 4 (1994) pp 200-206 Bloor, M I G and Wilson, M J ‘Complex PDE surface gcneration for analysis and manufacture’ Computing Vol 10 (199% pp 61-77
REFERENCES I
2 3
4
5
6
7
8
9 IO 11
152
Bloor, M I G and Wilson, M J ‘Generating blend surfaces using partial differential equations’ Comput. Aided. Des. Vol 21 No 3 (19891 pp 165-171 Hoschek, J and Lasser, D Fudamentais of Computer Aided Geometric Design A K Peters, Wellesley, MA (1993) Lowe, T W, Bloor, M I G and Wilson, M J, ‘The automatic design of hull surface geometries’ J. Ship Res. Vol 38 No 4 (1994) pp 319-328 Dekanski, C, Bloor, M I G and Wilson, M J ‘The representation of marine propeller blades using the PDE method’ .I. Ship Res. to appear (1995) Bloor, M I G and Wilson, M J, ‘Interactive design using partial differential equations’ in Sapidis, N (Ed.1 Designing Fair Curces and Surfaces SIAM, Philadelphia (19941 pp 231-25 1 Bloor, M I G and Wilson, M J ‘Using partial differential equations to generate free-form surfaces’ Cornput-Aided Des. Vol 22 (1990) pp 202-212 Vida, J, Martin, R R and Varady, T ‘A survey of blending methods that use parametric surfaces’ Cornput-Aided Des. Vol 26 No 5 (19941 pp 341-365 Brown, J M, Bloor, M I G, Bloor, M S and Wilson, M J, ‘Generation and modification of non-uniform B-spline surface approximations to PDE surfaces using the finite-element method, in Ravani, B (Ed.) Advances in Design Automation Vol 1 Computer Aided and Computational Design ASME (1990) pp 265-272 Prenter, Splines and VariationalMethods Wiley-Interscience, New York (1975) Burden, R L and Faires, J D Numerical Analysis (3rd edn) Prindle, Weber and Schmidt, Boston (19811 Ortega, J M and Poole, W G An Introduction to Numerical Methods for Differential Equations Academic Press, New York (1981)
Malcolm Bloor graduated with a BSc in mathematics from Manchester Unitersity in 1962. He obtained a PhD in applied mathematics from Manchester Uniwrsity in 1966. He was appointed as an assistant lecturer in the Department of Applied Mathematical Studies in 1964. He cuwently has a Chair in Mathematical Engineering at Leeds His research interests include computeraided design, fluid mechanics, industrial 1 mathematics, and biological mathemattcs
Michael Wilson graduated in 1980 with a BA in natural science from Cambridge University After research as pun of the Radio Astronomy Group in the Caoendish Laboratory, he obtained his PhD in I984 from Cambridge Unitersity. He joined the Department of Applied Mathematical Studies at Leeds University as a research fellow in 1983, and was appointed to a lecturing position in the same department in 1986. He is currently a senior lecturer. His research interests include computer-aided design, computational fluid dynamics, industrial mathematics and astrophysics.