COMPUTER AIDED DESIGN OF DEVELOPABLE SURFACES By Brian E. Konesky B.A.Sc. (Mechanical Engineering) University of British Columbia A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1993 © Brian E. Konesky, 1993 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Mechanical Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: Abstract The design of objects employing developable surfaces is of engineering importance because of the relative ease with which developable surfaces can be manufactured. The problem of designing developable surfaces is not new. Two space curves, defining the edges of the surface, are first created, then a set of rulings are constructed between the space curves under the constraint of developability. A problem with existing algorithms for designing developable surfaces is the tendency to introduce non developable portions of the surface; areas of regression. A more reliable solution to the problem of creating a developable surface is proposed. The key to the method is to define the developable surface in terms of a normal direc trix. The shape of the normal directrix defines the shape of the developable surface. Algorithms are defined to compute the shape of a normal directrix from a pair of space curves. A non-linear optimization technique was implemented to further refine the shape of the developable surface, but failed to yield satisfactory results. Other algorithms were also created to intersect adjacent developable surfaces and to generate the flat plate lay outs. The algorithms were implemented using the C++ programming language and the AutoCAD CAD package. Recommendations for further work are given. 11 Table of Contents Abstract List of Tables ix List of Figures x Nomenclature xiii Acknowledgement xvii 1 Introduction 1.1 1 Areas of Application Naval Architecture Aerospace Manufacturing Textiles 1.2 Definitions and Terminology 1.3 Methodologies 1.4 2 1.3.1 Kilgore’s Solution 1.3.2 Nolan and Clement Computer Approach 1.3.3 Normal Directrix Approach by Dunwoody . . New Approach Research Objectives Splines 14 111 3 2.1 Types of Splines Chosen 2.2 Uniform Parametric, Geometric and Non-Rational Continuity Requirements 18 2.3 Uniform Non-Rational Tension Catmull-Rom Spline 21 2.4 Uniform Non-Rational Beta-Spline 24 . 14 Creation of a Normal Directrix from Two Space Curves 30 3.1 Modern Approach Utilizing a Single Normal Directrix 31 3.1.1 Constraints Governing Modern Approach 32 3.1.2 Alignment of End Generators 33 Change in Angle With Respect to Unit Motion of a Control Vertex 34 3.1.3 The Alignment Process and the Concept of Mobility 36 Analysis of Results 38 Mobius Strip Demonstrating Flexibility 38 Controlling Alignment with Mobility 39 Problems arising when Surface has Small In and Out-of-plane Cur vature 3.2 41 Constraints Defining Modified Conventional Approach (Modified Nolan’s Approach) in Order to Create a Normal Directrix 42 3.2.1 Offset of Normal Directrix 47 3.2.2 Allowment of Different Number of Control Vertices 48 3.2.3 Present Problems with Current Solution and Improvisation Imple mented 3.3 49 Utilizing Modified Conventional Approach to Approximate a Single Direc trix 50 3.3.1 Equations Yielding Approximated Normal Directrix 52 3.3.2 Results and Present Problems 53 iv 4 Optimization of a Normal Directrix 60 4.1 Objective 60 4.2 Optimization Function 60 4.3 4.4 5 4.2.1 Integral Chosen to Minimize and Simplex Parameters Used 4.2.2 Present Problems 64 66 Downhill Simplex Method 67 4.3.1 Explanation of The Downhill Simplex Method 4.3.2 Results and Present Problems . . Examples 68 69 69 4.4.1 Example of Single Surface with Conical Properties 69 4.4.2 Example of UBC Series Demonstrating Present Problems 70 Intersection of Developable Surfaces and Flat Plate Layout 73 5.1 73 5.2 Intersection of Developable Surfaces 5.1.1 Derived Equations Yielding Intersections of Surfaces 74 5.1.2 Results and Present Problems 77 Conical Type Surface Tested 77 Criterion of Phantom Surfaces 78 Flat Plate Layout 5.2.1 5.2.2 6 . 80 Derived Equations giving Flat Plates Dimensions and Interplate Angles 80 Format of Output 81 Choice of Computer Language and CAD Program 83 61. 83 Computer Language Chosen 6.1.1 OOP 6.1.2 Portability to Different Platforms - Object Oriented Programming V 84 85 6.1.3 6.2 6.3 C++: Classes and Features 85 Class vector 85 Class matrix 90 Class Curve and Surface 91 Class ODE 93 Computer Platform Selected 94 6.2.1 PC 386/486 Environment 95 6.2.2 Computer Language Selection in this Platform 95 CAD Program Selected 96 6.3.1 CAD Program Environment and Open Architecture 96 6.3.2 Present Limitations 97 7 Demonstration Examples 8 98 7.1 Developable Mobius Strip 7.2 Simple Conical Developable 100 7.3 Arctic Fishing Vessel 100 7.4 UBC Series Fishing Vessel 100 98 Conclusions and Recommendations 105 8.1 105 Conclusions Normal Directrix from Two Space Curves Intersection of Developable Surfaces and Flat Plate Layout Implementation 8.2 105 . . 106 106 Recommendations 107 Bibliography 108 Appendices 110 vi A Mathematical Notation for Partial Differentiation 110 B Derivation of Developable Surface 111 B.1 Constraints which Define the Developable Surface 111 B.2 Proof Using Constraints 113 C Derivation of Rate of Rotation of Generator Differential Equation 115 D Tension Catmull-Rom Spline 122 D.l Phantom point, P(O), at 0 123 D.2 Phantom point, P(1), at n 125 E Beta-Spline 127 E.l Phantom point, P(0), at 0 128 E.2 Phantompoint,P(1),atn 131 F Derivation of Normal Directrix Control Vertices 134 G Modified Conventional Approach Derivation 139 H Relating the Space Curves to the Developable Surface 142 I Intersection of Developable Surfaces 146 J Derivation of Flat Plates 151 K O.D.E. Class, Adaptive Step Size and Runge-Kutta Method 153 L Non-Linear Optimization Downhill Simplex Method 155 M Codelisting 160 vi’ M.1 Class Tools Used . 160 M.1.1 vector.h 160 M.1.2 matrix.h 161 M.1.3 develop.h 162 M.1.4 ode.h 164 M.L5 vector.cpp 165 M.1.6 ma.trix.cpp 169 M.L7 solve.cpp 175 M.1.8 develop.cpp 179 M.1.9 ode.cpp 194 viii List of Tables 5.1 Flat plate data L.1 Values of variables used for Simplex Method . ix 82 157 List of Figures 1.1 Developable Surface Single Chine Hull 2 1.2 Manual method of generating developable surfaces 3 1.3 F117A stealth fighter 4 1.4 Definitions of Terminology of a Developable Surface 5 1.5 Kilgore’s Method of Creating a Developable Surface 6 1.6 Kilgore’s Manual Graphical Solution 7 1.7 Nolan’s Vector Description of Computer Solution 9 1.8 Normal Directrix Approach by Dunwoody 11 2.1 First Three Levels of Parametric Continuity 19 2.2 Comparing Geometric and Parametric Continuity 21 2.3 Catmull-Rom Splines 25 2.4 Beta Splines 29 3.1 Derivation of Developable Surface 33 3.2 Vector Locations and Corresponding Angles 34 3.3 Control Vertices and End Phantom Point 37 3.4 Robustness of Modern Approach Showing Mobius Strip 39 3.5 Mobius Strip 40 3.6 Modern Approach Showing Full Mobility of All Points 41 3.7 Provision Made for When Surface Very Near Flat 42 3.8 Relating Modified Conventional Approach Space Curves and Normal Di rectrix 43 x 3.9 In-plane and out-of-plane components 44 3.10 Positional Offset of Normal Directrix 48 3.11 Fanning Occurring When Developable Surface Nearly Flat 50 3.12 Conic-Like Section with Flatness Tolerance Specified at 0.001 53 3.13 Conic-Like Section With No Flatness Tolerance Specified 54 3.14 Developable Surface Procedure Comparison, tolerance 0.00 1 55 3.15 Developable Surface success, tolerance 0.01 56 3.16 Developable Surface Procedure Comparison, tolerance 0.01 57 3.17 Developable Surface success, tolerance 0.1 58 3.18 Developable Surface Procedure Comparison, tolerance 0.1 59 4.1 Relating the Distance Between Space Curves and Surface 4.2 Area calculated under space curves 65 4.3 Possible inherent failure due to certain geometry 67 4.4 Analogy of a Simplex For The Downhill Simplex Method 4.5 Conical-like developable surface 70 4.6 Developable Surface Failure 71 4.7 Developable Surface success, tolerance 0.001 5.1 Intersection of Three Developable Surfaces 5.2 Conical-type Surfaces Intersections • . . . 5.3 Conical-type Surfaces Intersections • • . • 5.4 Flat plate layout derivation 5.5 Developable surface and flat plate layout 81 6.1 Computer Platform Selected Initially 94 7.1 Developable mobius Strip 99 xi . . . 68 . 72 . 74 78 79 7.2 Conical-type Surfaces Intersections 101 7.3 Arctic Vessel Conventional Approach 102 7.4 Arctic Vessel Modern Approach 103 7.5 UBC series Vessel Conventional Approach 104 B.1 Derivation of Developable Surface 111 B.2 Derivation showing normal is invariant along a generator . . 112 C.l Vector locations and corresponding angles 115 G.l Orientation of Space Curves and Directrix 139 1.1 Intersection of Three Developable Surfaces 146 L.l Analogy of a Simplex for the Downhill Simplex Method 156 xii Nomenclature j,j 5 4 Dirac Delta Function. 0 Out of Plane Rotation of G’ With Respect to G. Change in Angle of Generator With Respect to Unit Motion of a Parameter De scribing Directrix. Rate of Change of Angle of Generator With Respect to Motion of Parameter De scribing Directrix. In Plane Rotation of G’ With Respect to G. a Parameter Describing Distance Along Generator. a 4 Incremental Position Along Generator. C: A Control Vertex that is Being Solved For. f(u) Position Along Left Adjacent Space Curve. Tangent Vector Along Left Adjacent Space Curve. Curvature Vector Along Left Adjacent Space Curve. g(v) Position Along Right Adjacent Space Curve. Tangent Vector Along Right Adjacent Space Curve. XIII Curvature Vector Along Right Adjacent Space Curve. . dg:ri The Developable Surface Generator Vector. Change of Position of Generator With Respect to Motion of Parameter Describing Directrix. d g 2 en Rate of Change of Generator With Respect to Motion of Parameter Describing dadt Directrix. dgen The Generator Slope Vector. Normal Vector For Left Adjacent Space Curve. Normal Vector For Right Adjacent Space Curve. n(t) Normal Vector From Surface ‘it Number of Control Vertices For Normal Directrix. n Number of Control Vertices For Left Adjacent Space Curve. n, Number of Control Vertices For Right Adjacent Space Curve. N Unit Normal Vector. p(t) Point On the Directrix. The Directrix Tangent Vector. - --- The Directrix Curvature Vector. s Scalar Value Along Generator xiv q(t, s) Point on Surface Determined by Parametric Values t and Scalar s t Independent Parametric Variable Describing Position Along Directrix. T Unit Directrix Tangent Vector. Change of Position of Unit Tangent With Respect To Motion of Parameter De scribing Directrix. u Dependent Parameter Describing Position Along Left Adjacent Space Curve. Incremental Position Along Left Adjacent Space Curve. v Dependent Parameter Describing Position Along Right Adjacent Space Curve. Incremental Position Along Right Adjacent Space Curve. xv A. Einstein Keep it simple: as simple as possible, but no simpler ENGINEERING The scientist analyzes what is. The engineer creates what has never been. The engineer scientist analyzes what is Imagines what should be Creates what has never been Analyzes the results of the creation. By Gunnar Schienius xvi Acknowledgement I would like to express my sincere thanks to Dr. A. B. Dunwoody for his enthusiastic support, supervision and direction on this project. I would also like to acknowledge the financial support from the NSERC of Canada whose financial support made this project possible. xvii Chapter 1 Introduction 1.1 Areas of Application Developable surfaces form a special class of surfaces which are very useful in many prac tical situations. Developable surfaces have many applications. A few applications are cited below. Naval Architecture Up until recently in history ship hulls were made of various types of wood, which could be easily worked into any desired shape. “The hull of continuous, homogeneous, testable sheet material is inherently stronger and lighter than the structure of small pieces of wood. If a skin of sheet material can be designed for low labour cost in construction, simple tools, and economy in repair, its engineering superiority and eventual economic advantage make it at once preferable to planks” [15]. To lower the labour cost in construction and repair even further, the skin of sheet material must be developable. At the same time, the hydrodynamic performance of the hull must be competitive with that of the best possible in compound hull construction [15]. Making the construction surface developable was therefore a desired criterion if possible. Figure 1.1 below shows a developable surface single chine hull of a fishing vessel. Figure 1.2 demonstrates how some manufacturers today create developable surface hulls for fishing vessels. 1 Chapter 1. Introduction 2 Figure 1.1: Developable Surface Single Chine Hull Aerospace Today’s modern aircraft are more complex in design and cost savings in manufacturing is crucial. Aircraft fuselages, wings and other smaller components can be produced by developable surfaces if considered in the design stage. One of the most recent “high-tech” declassified aircraft is the “wobblin-gobblin”, ie. the “Fl 1 7A”, (see figure 1.3) which has a low radar profile signature by using flat plates, is yet another example of using developable surfaces in an ingenious manner. Manufacturing In the area of manufacturing two new areas involving the developa bility criterion are in the application of peripheral milling and rolling. Each pass of an end mill cutting peripherally follows a developable surface. Chapter 1. Introduction 3 Figure 1.2: Manual method of generating developable surfaces Another area is in the application of rolling. Developability must also hold for this process. Both of these applications will not be discussed but only cited as other examples where developable surface criterion must hold due to the geometry of the application. Textiles Another example of developability is in the textile industry. Clothing is manufactured from fiat material and folded to conform to the appropriate geometry. Sails for sailing vessels is yet another application where developability must hold. 1.2 Definitions and Terminology The developable surface is a subset of the class of ruled surfaces. The definitions of a ruled surface and a developable surface are as follows: Ruled Surface: A ruled surface is defined as “The locus of a line, called a generator, Chapter 1. Introduction 4 Figure 1.3: F117A stealth fighter whose direction is determined by successive values of a parameter, moving continuously along a curve (a directrix) and intersecting that directrix at an angle other than zero.” [15]. Developable Surface: A developable surface,also defined by Kilgore and from Kreysig, is ‘A ruled surface having the same tangential plane on one and the same generator” [15]. Figure 1.4 below illustrates the terms introduced similar to the definitions presented by G.D. Aguilar[2j. The directrix must lie in the surface and that each plane tangential to the surface must also be tangential to the directrix. It must also be noted that with two space curves a ruled or compound surface may exist but no developable surface may be possible. Another theorem should also be noted that, “If two space curves lie in any developable surface, they lie in one and only one such surface” [15]. If the generators do not intersect anywhere, then the surface is developable. Chapter 1. Introduction Space Curves 5 Generator or Ruling Edge of Regression Figure 1.4: Definitions of Terminology of a Developable Surface In figure 1.4 the areas where the generators overlap are known as areas of regression. The boundary of an area of regression outlined in figure 1.4 is called the edge of regression [71. 1.3 Methodologies All existing methodologies for the design of developable surfaces start with the definition of two edges of the surface. Then, a set of generators is fit between those edges to define a developable surface. 1.3.1 Kilgore’s Solution The general method of matching developable surfaces to desired curves is an arduous task. The method assumes that the developable surface is either conical, cylindrical, or a combination of both. So, one tries a succession of surfaces until one is found to fit Chapter 1. Introduction 6 approximately. If the designer cannot find an exact solution his usual solution is to alter the original curve to fit the surfaces haphazardly [151. See figure 1.5. Kilgore’s Technique \‘ Figure 1.5: Kilgore’s Method of Creating a Developable Surface Kilgore examined this unique art and proposed a method for direct generation of developable surfaces from given beginnings. This method provides a manual graphical solution of surfaces to fit the curves, rather than to require alteration of the curves to fit the surfaces. Kilgore’s manual graphical solution is described in his paper[15] as well as a compre hensive description of the procedure is given in the Principles of Naval Architecture[1 11. A sample manual graphical procedure is displayed below which was extracted from the Principles of Naval Architecture[ll]. See Figure 1.6. One can see that manual graphical solutions are only as accurate as the skills and precision of the naval architect. This poses problems in error of the final solution as to its accuracy. C) 0 WI CD cm -a CD I-, oq r ‘ a p r pçn I I I -p ‘1 0 S n 0 C 11 I 4 3 0 0 I, .4 C -4 rn C n -I -I rn n x r > z 0 ‘1 c-fl ni I z -v n n I Chapter 1. Introduction 8 This led to implementing a mathematical solution once computers had evolved to the point where it was a viable alternative. The first computer solutions which had a moderate success were those by Nolan [17] and Clements [7]. 1.3.2 Nolan and Clement Computer Approach One of the first well known published works involving a computer-aided approach to developable surface design was by T. J. Nolan [17]. In his paper he emphasized that a mathematical approach utilizing a computer proved to yield a substanitial increase in speed and precision for calculating a developable surface. Nolan noted that an infinite number of surfaces can be found to span the curves but the developable surface is unique in that it requires the minimum strain energy of flexture and that in a developable surface bending is restricted to nonintersecting axes lying in the surface so that section moduli and bending stresses are minimized for any given radius of curvature. As a result, a developable surface can be formed elastically from a plane sheet, while the surface fitting the same pair of curves but having compound curvature must undergo a costly plastic forming process. Nolan defines a developable surface as, “a developable surface spanning a pair of curves in space may be defined as the locus of straight lines or “rulings” which represent the line of contact of a plane which is tangent to the surface. The rulings are neutral axes of bending and must not intersect within the surface” [17]. Nolan’s approach utilizes a Theilheimer spline interpolation and a vector represen tation to create a mathematical description for a developable surface. The approach is relatively simple, involving a representation of the tangents, normals and rulings in vec tor form. He iteratively solves his mathematical model to yield a zero angle for the cross product of the two normals of the space curves. See figure 1.7. Nolan’s vector approach calculates the normal of each space curve as N 1 and N 2 = = 1? x T , 1 R x T, where R is the ruling and T 1 is the tangent calculated at that point Chapter 1. Introduction 9 Ni df du nerator Cuef(u 9=( ; 9 f) xp(t) eg(v) Targent Plane Figure 1.7: Nolan’s Vector Description of Computer Solution along space curve 1 and T 2 is the tangent calculated at a point on space curve 2. He then uses one of the constraints of a developable surface, namely that N 1 xN 2 = 0. This approach is moderately successful in that for very simple surfaces it may yield a resulting developable surface. However, all too often the rulings either cross or “fan” yielding an unrealistic or nonusable surface. Also, the Theilheimer spline interpolation is not of parametric form resulting in singularities which may occur in a three dimensional representation of the surface. Clement’s Solution to the problem was published approximately ten years later uti lizing cubic spline functions, again in non-parametric form[7]. He states, “Between each pair of chine lines that ruled surface generated which has the same tangent plane at all Chapter 1. Introduction 10 points of each generator or ruling line. A procedure based on the multiconic development of a surface is used to modify the given chine lines to ensure that no ruling lines inter sect at a point within the surface. The result is a developable surface” [7]. In addition, Clement’s approach generated tables of offsets. This approach, like Nolan’s, also had problems arising at singularities and generators crossing in areas of regression on the surfaces. 1.3.3 Normal Directrix Approach by Dunwoody The approach taken by Dunwoody is unique in that it defines a developable surface by means of a normal directrix and an initial generator. The directrix is modelled by a parametric cubic spline; initially a uniform non-rational B-spline was used. Information about the spline in the parametric form is the position in space at a parametric value t, its tangent and its curvature. A start generator direction is also needed. The differential equation derived must retain the constraints which define a devel opable surface. Refering to figure 1.8 to show the geometry, the following constraints are shown and the resulting differential equation is formed. Details of the proof can be seen in appendix B of this thesis. From figure 1.8 the following vector definitions are in order: = generator vector (1.1) = derivative generator vector (1.2) = directrix tangent vector (1.3) = directrix curvature vector (1.4) = step increment (1.5) Chapter 1. Introduction 11 AT g g + /÷AT AT ‘h dp÷dpATT E+EAT dt 2 dt Figure 1.8: Normal Directrix Approach by Dunwoody Three constraints are necessary for a surface to be developable. They are stated as follows using the above nomenclature: 1. The normal directrix and generator vectors must be perpendicular. g dp Differentiation with respect to t yields g 2. The vector is of unit length. 77=’. Differentiation with respect to t yields Chapter 1. Introduction 12 -; g 3. The normal is invariant along a generator. ( X _; g) - = 0 Combining the constraints yields the following differential equation describing the next consecutive generator: (;-‘\ g_ dt — 16 1 d p dt dp dt Integrating this differential equation using such integration routines as Runge-Kutta fourth order forces a solution to be output. This approach will yield a particular solution. From this stage of the analysis, given a directrix and a starting generator of unit length, a unique developable surface results from the differential equation described above. There are limitations with this technique. It is not yet in a form which would prove to be of any practical use. Further constraining is necessary in order to control the behavoir of the developable surface. This approach is where the present research project started and expanded implement ing new terminology and concepts which will be presented in detail. 1.4 New Approach Research Objectives The new approach research objectives were based on expanding the work initiated by Dunwoody utilizing the normal directrix approach. This approach is unique in that it Chapter 1. Introduction 13 will always yield a solution. In this state, however, it is not very useful from a practical engineering view point since it does not match two space curves. This leads to the research objectives presented in this thesis to further refine this technique. • The first research objective was to develop an algorithm in order to find a normal directrix such that the resulting developable surface lay close to two space curves, representing desired edges of the developable surface. • Another research objective was to create an algorithm to intersect developable surfaces and to generate the flat plate layouts and angles. • The final objective was to implement these algorithms using a modern computer language and a popular CAD package in order to assess the practicality of the approach. Chapter 2 Splines This chapter is included because the material covered on splines contains the necessary background in order to derive the additional tools and equations for developable surfaces. The two types of splines discussed in this chapter were selected because of their particular characteristics which proved to be useful for a designer. 2.1 Types of Splines Chosen When considering using a mathematical representation for space curves one can classify them as of either parametric or non-parametric form. Non-parametric forms are used extensively in various fields of mathematics and engineering. Non- parametric curves can be further categorized as either explicit or implicit. The explicit non parametric form is usually expressed in the following form: y = f(x) where, x = independent variable f(x) = function of independent variable y = dependent variable 14 Chapter 2. Splirzes 15 In this form multiple-valued or closed curves cannot be expressed. To overcome this form, one usually uses the implicit form of the non parametric curve in the following typical form: f(x,y) —0 where, f(x,y) = A function of both x and y One typically calculates a point on the curve by calculating the roots of the equation. The approach can sometimes prove to be fairly computationally expensive. Implicit formulation is a very common form of non-parametric polynomials. Many formulations used in engineering require higher than third order thereby making computations even more expensive when solving for roots of the equation. The non-parametric implicit formulation presents difficulties when being applied in defining such three dimensional objects as ship hull curves and surfaces. One typical problem that arises is when a vertical slope is encountered along the curve or surface resulting in an infinite numerical value. An infinite number cannot be used in a numerical calculation when using a computer. Another problem arising in non parametric implicit form is that the positions are not distributed evenly along the curve or surface. This poses a problem when trying to present the curves or surfaces graphically on computers [19]. Parametric curves solve the problems presented above and are suitable for represent ing closed curves and curves with multiple values of an independent variable. Parametric curves are also axis independent. Parametric curves replace the use of geometric slopes Chapter 2. Splines 16 (which may be infinite) with parametric tangent and curvature vectors, which are never infinite. In parametric form a curve is usually represented by a piecewise polynomial. Each segment of the curve is given by three functions x(t), y(t), z(t), which are polyno mials in the parameter t. [12] For example: F(t) = [x(t),y(t),z(t)] After determining that parametric curves would be used in this work, the order and desired characteristics of the polynomials were investigated. Considering the types of applications and desired flexibility of the types and characteristics of curves desired, a fairly exhaustive investigation of various types of splines was conducted. Special cubic polynomials derived in the format pioneered by Barsky[3], DeRose[9], Forrest, Coons, Bezier and furthered by others were decided upon. The initial stages of development of the Basis functions to taylor the desired behaviour of the splines yield cubic polynomials. Cubic polynomials are most often used because lower-degree polynomials give too little flexibility in controlling the shape of the curve, and higher-degree polynomials can introduce unwanted “wiggles” and also require more computation. No lower-degree representation allows a curve segement to interpolate (pass through) two specified end points with specified derivatives at each endpoint. Given a cubic polynomial with its four coefficients, eg. f.’_a,43 x, — 1 j2 +VX +c+ four knowns are used to solve for the unknown coefficients. For example, the four knowns might be the two endpoints and the derivatives at the endpoints. Other knowns might be slopes or additional points[12]. It should also be noted that parametric cubics Chapter 2. Splines 17 are the lowest-degree curves that are nonplanar in 3-D (three dimensions). You can see this by recognizing that a second-order polynomial’s three coefficients can be completely specified by three points and that three points define a plane in which the polynomial lies. Higher-degree curves require more conditions to determine the coefficients and can “wiggle” back and forth in ways that are difficult to control. The parametric curves used in this thesis are given in terms of their degree n, which is fixed at 3. Much research has been done by such modern pioneers as Barsky[3j, who developed Basis functions, and appropriate nomenclature on the various levels and types of curve continuities. No detailed analyses of the derivation of the splines used in this thesis will be discussed in substantial detail since this work has already been done by such authors as those previously cited. Only enough explanation of the nomenclature used in this thesis to familiarize the reader with the concepts and characteristics of the various forms of the splines used will be discussed. Another point to mention is that local control of the 3-D curves was desired so that a curve segment is completely controlled by only four control vertices; therefore, a point on a curve segment can be regarded as a weighted average of these four control vertices. The parametric splines used in this thesis are preseilted in the following form: 1 a a 3 a 4 2 a 1 b P(i + t) = 3 2 b 3 b 4 b 1 F_ F, (2.1) i 1 C 2 C 1 d d 2 3 C 4 C 3 d d 4 Chapter 2. Splines 18 where, P(t) = [x(t), y(t), z(t)j = P(t) (2.2) and Point Tangent Curvature 2.2 T(t) = C(t) [t3 = = = &t) = = 2 t i] t [3t2 [] [ 2t 1 0] [6t 2 0 0] [ } [] [1 [1 (2.3) (2.4) (2.5) Uniform Parametric, Geometric and Non-Rational Continuity Require ments One of the important properties discussed in such fields as finite elements and computer aided geometric design is of the mathematical techniques of shape representation. It is termed continuity. Continuity can be described as the highest level of differentiation which is continuous [3]. The types of continuity can be further categorized. Four types of continuity are considered in this analysis. Each type of continuity will be explained briefly. The types of continuity chosen can be expanded but only two of what was thought to be generally the most useful were selected at this stage of the research. The first type of continuity requirement considered is whether or not the splines used in the analysis are either uniform or non-uniform. Since the splines are parameterized, and uniform parametric splines can be expressed in a pseudo standard format, these types of splines appeared to be a likely choice. Non- uniform splines are not able to be expressed in the format that was adopted in this reseach at this time. Another type of continuity requirement which was desired was parametric continu ity [18]. Parametric continuity can be explain quite briefly with the aid of Figure 2.1. If t1 derivative vector of two cubic curve segments are equal (ie. their direction and the n Chapter 2. Splines 19 magnitudes are equal) at the segments’ join point, the curve has nthdegree continuity in the parameter t, or parametric continuity[12]. One would then state that if the direction and magnitude of [P(t)] through the th derivative are equal at the join point, the curve is called C’ continuous. Figure 2.1 shows a curve segment S joined to three differ ent curves with three different degrees of continuity ascertained by the superscript above the C. One should note that a parametric curve segment is itself everywhere continuous; the continuity of concern here is at the join points[12j. y(t) 2 IC x(t) Figure 2.1: First Three Levels of Parametric Continuity If two curve segments join together, the curve has GO geometric continuity. If the directions (but not necessarily the magnitudes) of the two segments’ tangent vectors are equal at the joint point, point the curve has G’ geometric continuity. In computer- aided design of objects, G 1 continuity between curve segments is often required. G’ continuity means that the geometric slopes of the segments are equal at the join point. For two tangent vectors TV 2 to have the same direction, it is necessary that one be 1 and TV a scalar multiple of the other. We then state the relationship that TV 1 k> 0 [3j[12j = k TV 2 with . Chapter 2. Splines 20 One should note that in general, C’ continuity implies G’, but the converse is gen erally not true. G’ continuity is generally less restrictive than is C , so curves can be 1 1 but not necessarily C G . However, visually, join points with C 1 1 continuity will appear just as smooth as those with C’ continuity as can be seen in Figure 2.2.[12j. 2 and are identical except Q,, Q, Q join at the point P for their tangent vectors at P . Q, and Q2 have equal tangent vectors and hence are 2 both G’ and C’ continuous at P . Q, and Q have tangent vectors in the same direction 2 but Q has twice the magnitude so they are only G’ continuous at P [12] 2 In Figure 2.2 curve segments , , Another type of continuity which one may desire is Rational Continuity [13]. Rational continuity can be simply defined for “general rational cubic curve segments as ratios of polynomials: — — X(t) 147(t)’ — — Y(t) z( 147(t)’ — — Z(t) 147(t) 2 6) . where X(t), Y(t), Z(t) and W(t) are all cubic polynomial curves whose control points are defined in homogeneous coordinates. We can also think of the curve as existing in homogenous space as: Q(t) = [X(t)Y(t)Z(t)TV(t)] (2.7) As always, moving from homogeneous space to 3 space involves dividing by W(t). Any non rational curve can be transformed to a rational curve by adding W(t) = 1 as a fourth element” [12]. No splines were used in this thesis at this stage which included Rational continuity. In future work this type of continuity requirement may be implemented if requested. For further reading one may refer to either Barsky and Hohmebar [13] or Foley [12j Chapter 2. Splines 21 y(t) ç,P Figure 2.2: Comparing Geometric and Parametric Continuity 2.3 Uniform Non-Rational Tension Catmull-Rom Spline The uniform non-rational Tension Catmull-Rom spline was chosen because it exhibits severa’ useful features the designer may require [8] [91. First, it is an interpolating spline, meaning that the curve passes through the points (control vertices). Second, it is in parametric form, meaning that it does not encounter singularities, only variations of vector magnitudes. Third, it exhibits desired parametric and geometric continuity requirements. Fourth, it has a global tension parameter which can further control the shape of the desired curve. The uniform non-rational Tension Catmull-Rom spline is easiest described in vector matrix form. In this form it is a relatively simple task to imbed into a program. The vector-matrix format is in a form in which not only the position along the curve can Chapter 2. Splines 22 be calculated but also the tangent, curvature and other vector specific relations desired. This vector-matrix format is now almost a standard form in which these types of splines are presented. This pseudo-standard form shows that the spline has a local influence of the control vertices as any chosen position. At any one particular location along the curve the control vertices are only influenced by the previous one and the next two. This is shown in the “[P]” vector of the vector-matrix form. The “[P]” vector shows that at a particular position along the curve the only infuence is from P_ , P, P:+iandPj+ 1 . 2 Taking into account the end conditions of the splines also had to be considered. The approach taken was to create Phantom points. Given the control vertices vector below, Phantom end conditions were formulated. Detailed analysis of the derivation can be referred to in the appendices D.1 D.2. The vectors located below show the indexing of the control vertices and how end conditions are treated (phantom points). 1 PiPi 1 Pi+ 2 Pi+ For the initial condition P 0 the control vertices vector has the following form D.l: Chapter 2. Splines 23 Pi Pi 1 Pi+ i+2 For the end condition P,_ 1 the control vertices vector is in the following form D.2: Pi Pi 1 Fi+ 1 Pi+ The following page shows the vector-matrix form of the uniform non-rational Tension Catmull-Rom spline. The spline is parameterized with respect to the parametric variable t and has a global tension variable, 6. —2.0/3 4.0 1 P(t) = 2 t’ 3 t t — 4.0/3 2.0/3 2.0/3 2.0/3 — 6.0 6.0 — — 4.0 2.0/3 4.0/3 —2.0/3 1 —2.0/3 0.0 2.0/3 0.0 0.0 2.0 0.0 0.0 Chapter 2. Splines 24 1 PiPi 1 Pi+ 2 Fi+ = Tension parameter To give the reader a good comparison as to the behaviour of the Tension Catmull Rom spline the following figures are included. The first figure, figure 2.3(a) shows the spline with a tension value of 0.5 shown relative to interconnected line segments. For the tension Catmull-Rom spline with a value of 0.5 this defaults to a traditional cardinal spline. The next figure, Figure 2.3(b) shows little difference but is relaxing the spline tension when given a tension value of 1.0. Finally the Tension Catmull-Rom spline in Figure 2.3(c) shows how it nearly contours to the inter-connected line segments shown in the figure. If the tension parameter is given a value of 0.0 then it becomes line segments. 2.4 Uniform Non-Rational Beta-Spline The uniform non-rational Beta-spline was chosen to provide other useful features. First, it is an approximating spline meaning that the curve passes near, not through, the control vertices. It also exhibits the convex-hull property which is also shared by the B-spline [18] and the Bezier [4] spline. Second, the Beta-spline is derived parametrically so it too does Chapter 2. Splines 25 (a) Tension at 0.5 Catmull-Rom (b) Tension at 1.0 Catmull-Rom (c) Tension at 0.1 Catmull-Rom Figure 2.3: Catmull-Rom Splines Chapter 2. Splines 26 not have any singularities occur. Third, it also retains desired parametric and geometric continuity requirements. And, fourth, it has global bias and tension parameters which further enable the designer to better adjust the spline [3]. The vector-matrix form of the Beta-spline is shown on the following page: P(t) 1 2 t —2.0/3? 2.0(132 + i3? + /3? + j3) 6.0/3? —3.0(132 + 2.0/3? + 2.0/3?) —2.0(132+/3? +/31 + 1.0) 2.0 3.0(132 + 2.0/3?) 0.0 6.0/3k 0.0 + 4.0,3? + 4.0/3i) 2.0 0.0 6.0(/3? (132 1] /3k) —6.0/3? 2.0/3? t’ — 1 PiPi 1 Pi+ = i3 + 2.0/3? + 4.0/3? + 4.O/3 + 2.0 /32 = Bias = Tension Chapter 2. Splines 27 Like the Tension Catmull-Rom spline the Beta-spline exhibits the same localized control vertices influence. Again, if one were to choose a particular position along the curve the closest control vertex would only be influenced by the preceding one and the next two, eg. [P_ , P, P+i, P÷ 1 1. 2 Taking into account the end conditions of the spline also had to be considered in the same fashion as the Catmull-Rom spline. The same approach was taken to create the Phantom points. Given the control vertices vector below, Phantom end conditions were formulated. Detailed analysis of the derivation can be referred to in the appen dices E.1 E.2. The initial condition control vertices vector for the Beta-spline at P(O) is: (P2+4.O3?+4.OI31 )P +2.oP+i 1 6—2.Of3 Fi 1 Pi+ The end condition control vertices vector for the Beta-spline at P(n-l) is: 1 PiPi 1 Pi+ 2 +(4.0+4.0/31+/32 )P 2.013?P -i-i t 3—2.0 Chapter 2. Splines 28 The first of the Beta-spline figures shown below reveals the spline compared to inter connected line segments. The first figure, figure 2.4(a) shows the curve near the control vertices, exhibiting the characteristics of an approximating spline. The bias is at 1.0 and the tension is at 0.0. With these values the Beta-spline degenerates to a B-spline, which exhibits first and second order parametric continuity. The next figure, figure 2.4(b), shows the Beta-spline with a Bias of 1.0 and a tension at 25.0. At these values one can see how this spline can also resemble inter-connected line segments if the tension value is increased substantially more. In figure 2.4(c) one can see how the curve behaves if the tension parameter is given a negative value of about -0.05. Finally, if the bias is changed, as in figure 2.4(d), to a value of 1.5 and the Tension is left at a value such as 0.0 the curve exhibits the following behaviour These figures shown above try to give the reader some idea of the capabilities of these type of splines and how one can use the features each spline possesses. These are but a few of the types of splines which are now being developed. Each of these types of splines has features which the reader should be aware of in order to maximize the benefits they have to offer. Chapter 2. Splines 29 (a) Bias 1.0 Tension 0.0 Beta (b) Bias 1.0 Tension 25.0 Beta (c) Bias 1.0 Tension -0.05 Beta (d) Bias 1.5 Tension 0.0 Beta Figure 2.4: Beta Splines Chapter 3 Creation of a Normal Directrix from Two Space Curves The normal directrix approach will always yield a smooth developable surface in the vicinity of the normal directrix, so long as the normal directrix is itself smooth. Unfor tunately, it is not always clear what shape the normal directrix should take in order that the resulting surface should meet the requirements of the designer. With respect to the requirements of the designer, the definition of a developable surface from two edges is superior. The objective of the present work is to create an algorithm which will shape a developable surface to lie close to a pair of edge curves. The developable surface will be specified in terms of a normal directrix in order to ensure that the surface is smooth. It is not expected that the developable surface will contain the two edge curves, only that it will be close to the two curves. The creation of a normal directrix from two space curves follows from these criterion: 1. A normal directrix can be computed for any developable surface by starting at one point on the first generator, then constructing a curve which lies within the surface and is perpendicular to all generators. In addition, extra construction tools, in the form of differential equations, were created in order to better control the normal directrix solution. 2. Once a normal directrix has been computed, it can presumably be smoothed out to yield a smoother developable surface without departing greatly from the original equation. 30 Chapter 3. Creation of a Normal Directrix from Two Space Curves 31 3. Nolan’s approach of matching the cross products between the generator and the tangents to each of the edge curves can be expressed in terms of a differential equation. 4. The curve of the normal directrix can also be expressed as a differential equation, to be solved in conjunction with the differential equation for the set of generators. 5. If the normal directrix is to be described by a spline with n control points, then the values for those n points can be computed to yield a spline which lies close to the normal directrix derived from the differential equations. The approach taken to try to match a developable surface to two edge curves re sulted in formulae modelled by differential equations. The differential equation version of the conventional method, named the modified conventional approach, created by Dunwoody and Konesky was used as an initial guess in order to utilize additional differential equations to solve for directrix control vertices. 3.1 Modern Approach Utilizing a Single Normal Directrix A normal directrix can be computed for any developable surface by starting at one point on the first generator, then constructing a curve which lies within the surface and is perpendicular to all generators. This very powerful technique was created by Dunwoody which is termed in this thesis as the modern approach. In addition, extra construction tools, in the form of differential equations, were created in order to better control the normal directrix solution. This modern approach, given a normal directrix and a start generator position, will always force a developable surface to be created. This solution is underconstrained, however, and further refinement was deemed necessary in order to better control the Chapter 3. Creation of a Normal Directrix from Two Space Curves 32 behaviour of the function. The formulation is in vector differential equation form and uses the ODE class integra tioll routine which implements Runge-Kutta 4th order. Constraints for this developable surface differential equation are given in the next subsection. For a more thorough anal ysis one can refer to Appendix B for the derivation. Constraints Governing Modern Approach 3.1.1 The differential equation developed by Dunwoody [10] termed the modern approach in volves three constraints which define a developable surface. They are as follows: 1. The generator must be of unit length ie. 2. The vector normal is invariant along a generator ie. (gx).=0 3. Vectors must be perpendicular ie. dg dt dp_ dt d p 2 g 2 dt The three constraints can be seen in Figure 3.1 below in vector form. Using the three constraints which define a developable surface yields the differential equation in simplified form: dg dt (. — — 2 dt \dt g,“dP dt) Chapter 3. Creation of a Normal Directrix from Two Space Curves 33 g dt ‘1Illi: dt 2 dt Figure 3.1: Derivation of Developable Surface This equation works fine as is but further constraining is necessary in order to make this solution practical. For example, given the theory presented so far we can generate a developable surface along a normal directrix given an initial generator position. However, more realistically, one would also want the surface to end up in alignment with a desired final generator position. We now move to the next stage in the development of alignment of the end generators. 3.1.2 Alignment of End Generators Once a normal directrix has been computed, it can presumably be smoothed out to yield a smoother developable surface without departing greatly from the original equation. In addition, a much more realistic and desireable condition is where the user gives a normal directrix, a start generator position and a final generator position and the configuration adjusts itself to conform to the newly added constraints. Chapter 3. Creation of a Normal Directrix from Two Space Curves 34 Change in Angle With Respect to Unit Motion of a Control Vertex The problem was approached by looking at the problem using the perspective of observing the location of the directrix, the position of the generator and separate the vectors into inplane and out-of-plane components. Another differential equation was formulated which accumulated information of the rate-of-rotation of a generator with respect to motion of the control vertices along the normal directrix . This formulation, as will be shown, proves to be very useful in storing information about the surface and how to correct accordingly to match the end generator. Figure 3.2 shows the relation of the original directrix and the corrected directrix as well as the amount of change that is necessary. N’ T — ___. ———Original Directrix Modified Directrix Figure 3.2: Vector Locations and Corresponding Angles A detailed analysis of the derivation of the differential equation can be referenced in Appendix C if the reader wishes to look further. A brief summary of the highlights of the derivation is needed here in order to familiarize the reader with new concepts which are being introduced with the theory and the nomenclature of the user controllable design variables. The angle Theta, herein referred to as 9, is the out-of-plane rotation of G’ with respect to G. The angle Phi, is herin referred to as 4’, and is the in-plane rotation of G’ Chapter 3. Creation of a Normal Directrix from Two Space Curves 35 with respect to G. The variable “a” is a parameter describing a directrix control vertex component. The rate of rotation of a generator with respect to changes in one of the parameters of the directrix curve is written as: dadt The desired expression is the rotation of a generator with respect to changes in one of the parameters of the directrix, which is written as: dO da The desired expression can be defined in terms we can derive, namely: = N = -•N (3.2) where, N is the unit normal defined as: (3.3) Txg (3.4) but, ddO dtda rewriting gives, d 8 2 dadt (3.5) — — d(dg N dt”da 36 (3.7) — — d g 2 dadt dt dN dt 38 Using numerous identities and proofs, the user can look at the derivation in detail in Appendix C, the resulting differential equation which contains both the in and out-of plalle components result in the following formulation: — dadt — (“) 39 ) . Chapter 3. Creation of a Normal Directrix from Two Space Curves 36 Using several identities and constraints that have already been presented the final form simplifies down to the following relation: O(x) (2 dadt \dt dt’ &p (3.10) — The Alignment Process and the Concept of Mobility The alignment process involves a summation of the term from t 0.0 to t = N-i. The summation can be written in equation form as follows: dadt {6—1.O} N—1 1 Jo !dt 4 dadt = N—1 1 Jo 2 ( \dt X p 2 d dadt{5} ) dt (3 11 dt’ dadt{5+1.O} dadt {6+2.O} where dadt is defined as: a bcd efgh = (3.12) 2 2.Ot 1.0 0.0 3.0t i j m n k 1 6+1.0 o w 6+2.0 Chapter 3. Creation of a Normal Directrix from Two Space Curves 37 The alignment process may take several passes, ie. from 0 to N-i, a correction, then from 0 to N-i, etc. The procedure will be described and the successive passes usually reduce the error by one decimal place per pass (ie. iteration). After the first accumulation of information along the spline from t = 0 to t = N 1, a few constraints are added. The first desired constraint of the movement of the control vertices is that the end control vertices can only move in the direction along each corresponding end generator. The second from the end control vertices are constrained to move both in the direction of the end generators and in the direction of the start and end tangent vectors respectively. Both ends are calculated in the same way, so for a better understanding only the t=0.0 and t=l.0 end conditions will be explained. At t = 0.0 the control vertex located here only has the freedom to move in the direction of the starting generator. The end condition constraints are also influenced by the phantom end conditions in addition to depending upon the type of spline being used. For this analysis the degenerated version of the Beta spline to a B-spline will be used. I Pantom Point Figure 3.3: Control Vertices and End Phantom Point At the end of the spline the end two control vertices of the summed terms have to Chapter 3. Creation of a Normal Directrix from Two Space Curves 38 be modified because of the phantom end condition constraints. The end control vertices are influenced by the following condition which describes only the B-spline criterion: 2 2F — . From this relation we have to add 1 P of information and subtract again to itself for the P 2 control vertex of the 2 P + ’ component from itself. This is in essence the technique which is applied to the end conditions of the control vertices for each type of spline being used. The concept of mobility is quite simple and provides a further constraining on the behaviour of the alignment process. Each control vertex of are assigned a mobility value, which is a constant. A mobility value of 1.0 means that the corresponding vertex has full mobility and that it is free to adjust that control vertex as governed by the equation. At the other extreme a mobility value of 0.0 indicates that the vertex is not to be moved during each iteration of the alignment process. 3.1.3 Analysis of Results Many tests were done with the modern approach and then it was incorporated as a foundation to build upon. The approach was found to be very robust but specifically controlling its behaviour became the dominant problem. It was also noted to have one major instability. This occurred when the plate was very close to being perfectly flat. This was not deemed to be a very major problem since if the plate was flat, a solution existed already, thereby, no need of the modern approach would be required. Mobius Strip Demonstrating Flexibility As an example which demonstrates both the flexibility of the modern approach ie. the control vertices were not permanently fixed in position, (NB. only the end conditions were given a hard constraint). One challenge was to create a developable mobius strip. The mobius strip is a geometric anomally in that its edges are infinitely continuous, ie. Chapter 3. Creation of a Normal Directrix from Two Space Curves 39 if one was to follow an edge it would continue infinitely along the surface travelling along both sides. For clarity, the reader should refer to Figure 3.4 in this chapter. Figure 3.4: Robustness of Modern Approach Showing Mobius Strip In order to create the developable mobius strip the very end control vertices were fixed. As mentioned earlier a new concept was introducted. This concept was termed mobility and is explained in ‘the next section below with figures included to help clarify the explanation. Controlling Alignment with Mobility The concept of mobility was developed as a first tool to control the behaviour of the modern approach. It can simply be defined as local weighting of the change in angle with respect to unit motion of each control vertex. Each control vertex has information recorded about it by the function described in equation 3.10. When a mobility value of 1.0, unity, is assigned the class entity which retains the information about a particular vertex has 100% freedom to reposition the location of that particular vertex. In contrast, a mobility of 0.0, indicates that the vertex is not permitted to be moved at all. We can use a range of values for each vertex ranging from 0.0 to 1.0 if we wish to “fine tune” or Chapter 3. Creation of a Normal Directrix from Two Space Curves 40 more accurately control the behaviour of the iterative solving procedure in order to align with the end generator. The example of the developable mobius strip shown previously and now shown in two views in figure 3.5 was constructed with 6 control vertices. Each control vertex was assigned a mobility. The first and last vertex were given a mobility of 0.0, ie. no movement allowed, and the other four were given full mobility of 1.0, ie. full freedom to be corrected. (a) Developable mobius strip view 1 (b) Developable mobius strip view 2 Figure 3.5: Mobius Strip The developable mobius strip was one of the first examples where mobility was neces sary to create a specific type of shape. The same 6 control vertices used to construct the mobius strip were all initially given mobility of 1.0, and run to see what type of solution would result. The resulting figure 3.6, demonstrates what type of solution results when Chapter 3. Creation of a Normal Directrix from Two Space Curves 41 no constraining is implemented. Figure 3.6: Modern Approach Showing Full Mobility of All Points One should note that figure 3.6 is a valid solution. The end generators do align with the same line that passes through both the start and end generator positions. Clearly one can see that if “hard” constraints are not given, more than one solution can result. This provided us with insight in furthering the design analysis in that more than one solution could result; an entire family of solutions is possible with the single directrix approach unless further constraining is included. Problems arising when Surface has Small In and Out-of-plane Curvature In the modern approach a problem arises when the surface has small in and out-of-plane curvature , ie. the surface is almost flat. This problem is located where the alignment process takes place. When the accumulated values of the change in angle with respect to motion of the control vertices is very near zero, a subsequent calculation involves Chapter 3. Creation of a Normal Directrix from Two Space Curves 42 dividing the present angle that the end generator makes with the desired end generator. This results in division of a floating point value with one which is very near to zero. If the plate is very close to being flat the computer being used to do the analysis will yield a floating point error due to division by zero. A tolerance or threshold value was implemented which would make the resulting floating point operation equal to zero if the total in and out-of-plane curvature was less than 1.0 x 10_b. This was deemed necessary in order to improve the robustness of the algorithm. This results in the modern approach to accomodate when the surface is very near flat when initialized. A sample is shown in figure 3.7 where the surface is flat and a corresponding correct solution results. H Figure 3.7: Provision Made for When Surface Very Near Flat 3.2 Constraints Defining Modified Conventional Approach (Modified Nolan’s Approach) in Order to Create a Normal Directrix Nolan’s approach of matching the cross products between the generator and the tan gents to each of the edge curves can be expressed in terms of a differential equation. Nolan’s approach can give good approximations of developable surfaces. This approach expressed in differential equation form and with additional constraints could be used as an approximation for creating a close fitting normal directrix. Before a normal directrix can be computed, it must have information as to where Chapter 3. Creation of a Normal Directrix from Two Space Curves 43 it should be located relative to two space curves that contain a developable surface. Information in the form of modelling the Conventional Approach by differential equations lead to the formulation of Normal Directrix Control Vertices. This material is presented with the Modern Approach Utilizing a Single Normal Direct rix because it is used to derive the normal directrix control vertices. In figure 3.8 the two outer space curves, referred to as f(u) and g(v), are used in defining the modified conventional approach. In this approach the first two of the three constraints stated by Nolan [17] are exactly the same. The third definition is also used but is modified to include an additional term is formulated in order to determine where along the generator the normal directrix should lie. These steps are now explained in more detail. Ni Cuef(u) Cueg(v) Tangent Plane Figure 3.8: Relating Modified Conventional Approach Space Curves and Normal Direc trix In figure 3.8 the first two constraints are graphically shown stating that each space Chapter 3. Creation of a Normal Directrix from Two Space Curves 44 curve must have each tangent vector perpendicular to its normal. Also, the cross-products between the generator and the tangents to each space curve must be parallel to each other. We can express the first two constraints for each space curve as follows: where t (3.13) —Xgen = du (g(v)—f(u))=(g—f) = dg dv and ii;; — X (3.14) gen (3.15) Out of Figure 3.9: In-plane and out-of-plane components The third relationship is derived from the out-of-plane components and the in-plane components leading to their next respective locations along the space curves. The third relationship is stated in equation 3.16 and is explained here below. d2 — dt (genxni) du dg — — (gen — X ) 2 n . dt (3.16) Chapter 3. Creation of a Normal Directrix from Two Space Curves 45 The numerator dot product relation in equation 3.16 represents the out-of-plane com ponent. The denominator relation represents the in-plane component. Only magnitudes are needed since direction is constrained to being along each space curve’s tangent and normal vectors. At a particular position along the space curve the out-of-plane direction is located along the normal vector and curvature vector. For example, on the space curve f(u), i and are in the directions of out-of-plane curvature. For in-plane curvature one can relate the resulting cross product of the generator, §, with the normal, i, to get a vector in the appropriate in-plane direction. The tangent vector of curve f(u), and the vector, lie in-plane. From these relations a ratio of just the magnitudes would be simplest to equate with respect to the next parametric calculation along the space curve since the directions are constraints in the formulations. These relations are then equated to two space curves. By equating these relations to two space curves a known developable set of generators (or rulings) can be first approximated. The equat ing of the ratio of out-of-plane curvature to in-plane direction for two space curves f(u) and g(v) in equaton 3.16 is written as the definition just described. One should note that equation 3.16 is close to the equivalent approach described by Nolan. We now take into consideration a different number of control vertices for all of the space curves,and an approximated position of where the normal directrix should lie. We now show the key equations which yield the final result in the modified conven tional approach. For a more comprehensive derivation of this refer to appendix G. Initially we need to set up the three space curve relationships: p(t) = (1 — a)f(u) + ag(v) From the two space curves, f(u) and g(v), the generator is found as follows: (3.17) Chapter 3. Creation of a Normal Directrix from Two Space Curves (3.18) (g—f) = 46 Equation 3.17 is then differentiated with respect to t in order to establish a differential equation for the normal directrix intersection with the generator. Differentiating yields: dp dfdu dgdv da = Using the definition that, . = (3.19) 0 and using this on equatioll 3.19 results in the following: dp df du dg dv (3.20) = The normals of each space curve is also calculated using the following relationships: = df = x gen (3.21) (3.22) Relating the out-of-plane curvature with the in-plane curvature of the two space curves as mentioned previously we get: : _; (gen x ni). = . (3.23) (gen x n ). 2 Including the relationship of non-equal number of control vertices for the space curves: Chapter 3. Creation of a Normal Directrix from Two Space Curves dv dt 2n — nt du 2ri,,dt (1 — Rearranging these relationships to solve for , and 47 3 24 ij gives the following relation: —1 2f d du dt (genxni).— — — 2n + (.325 ) 2n dv2 (genx n ).— 2 dv —- dv dt — — = 2n ndu ndt t 326 dg —s dv du (l—a)—.gnj+ a--.gen.. —.--- gen gen 3.2.1 (3.27) Offset of Normal Directrix One feature which was found useful for the designer to have as a variable to adjust is the initial positional offset of the normal directrix. This relation is shown again here as equation 3.28: p(t) = (1 — a)f(u) + ag(v) (3.28) By allowing the initial position offset, a, of the normal directrix to be adjusted the user may be able to further fine tune the surface to desired specifications. If this is not of any concern it defaults to the parametric value, 0.5, which is the midpoint between the two space curves. Chapter 3. Creation of a Normal Directrix from Two Space Curves 48 f(u) g(v) Figure 3.10: Positional Offset of Normal Directrix 3.2.2 Allowment of Different Number of Control Vertices During the initial design stages one criterion was noted, namely that the number of control vertices for each space curve and the normal directrix may be desired to be different from each other. This was noted and a further relationship was established which allows for different control vertices for any of the curves. The relationship is presented as follows: = The number of control vertices of the normal directrix. = The number of control vertices of one space curve. = The number of control vertices of another space curve. Given the variables shown above and the two differential parametric index parameters, and , a very simple equation can be stated as follows: Chapter 3. Creation of a Normal Directrix from Two Space Curves 1.0 dv 1.0 du = 49 1.0 (3.29) — Equation 3.29 states that each parametric differental variable, and , contributes one-half times its total number of control vertices. Ignoring the total number of control vertices for every curve results in the equation being written as: 1.0 + (3.30) The resulting relationship for the space curves is as follows: dt = flj ndt (331) Further reading detailing the derivation can be referred to in G. 3.2.3 Present Problems with Current Solution and Improvisation Imple mented Tests including equation 3.16 proved promising. One problem noted, is shown in Fig ure 3.11 below when the surface is close to becoming flat. A phenomenon which we termed “fanning” occurs where the out-of- plane component becomes very small relative to the in-plane-component yielding undesirable characteristics. This is shown below in Figure 3.11. One can see from Figure 3.11 that intuitively if a plate is close to becoming flat, the generators should lie nearly perpendicular to the tangent vectors of each space curve, resulting in square flat plates for a flat surface. This problem was noted and ad dressed by incorporating an out-of-plane tolerance which would default to parallel lines when the surface fell below the user-specified tolerance. This is a temporary solution and further research is being directed to address this problem in future work. Chapter 3. Creation of a Normal Directrix from Two Space Curves 50 Figure 3.11: Fanning Occurring When Developable Surface Nearly Flat For a more detailed example showing where this phenomenon occurs in practice see the examples in Chapter 7, entitled “Demonstration Examples”. The three examples cited were a simple conical developable, an arctic series fishing vessel, and the UBC series fishing vessel. 3.3 Utilizing Modified Conventional Approach to Approximate a Single Di rectrix If the normal directrix is to be described by a spline with n control points, then the values for those n points can be computed to yield a spline which lies close to the normal directrix derived from the differential equations. From the equations 3.25, 3.26 and 3.27, the intersection location along the generators can be calculated. The desired number of control vertices for the normal directrix is still to be determined and may vary. The number of desired normal directrix control vertices is specified by the variable From the above mentioned relationships another differential equation is needed in order to solve for the location of the normal directrix control vertices. The relationship was found to be in the form of finding the minimum for an equation relating the normal directrix in terms of an error. This relation relates p(t), the known spline, and r(t), the spline control vertices we wish to solve for, is as follows: error N—i = Ip(t)—r(t)I d 2 t (3.32) Chapter 3. Creation of a Normal Directrix from Two Space Curves 51 We then differentiate equation 3.32 with respect to the control vertices and solve for the relation to determine a minimum. This is shown as follows: Oerror JN1 = 2 3- (p(t) — r(t)) dl (3.33) N—2 = YJ = 0 ‘(p(j+s)-r(j+s))dS 0 (3.34) . (3.35) Details of the analysis can be seen in Appendix F. The initial relation showing the equation is as follows: T 1 Cj_ ij—1 8 s2 [A] = ij+1 6 S ij+2 8 1 [ . 1 ] [A] 1 cj+ 2 + 3 C 1 Cj_ ij—1 6 2 C Cj+l 2 C, dS (3.36) [ —1 N—2 = (ii)] N—2 1 r(j+s) [3 2 s 1 ] ‘ i 5 i [A] d3.37) ij+1 5 ij+2 5 Chapter 3. Creation of a Normal Directrix from Two Space Curves 3.3.1 52 Equations Yielding Approximated Normal Directrix Collecting equations 3.25, 3.26, 3.27, equation 3.37 and presenting them below shows how all of the variables are sequentially coupled and integrated using the class ODE (Ordinary Differential Equation) solver. —1 ______ 2 d du dt (genxn ) 1 .— — — + 2n 3 38 ) 2n . 2 dv (genxn ) 2 .— dv dv dt — — 2-i,, nt ndu ndt 339 - cia d1 — (1 — g du dv a)— gen + a— gen (3 40) 1 Cj_ sij—1 C, N—2 = 1 Cj+ [(ii)] —1 N—2 yj 1 r(j+s) [3 2 [A] d3.41) sij+1 ij+2 5 4 In the computer program all of the terms are collected as the ODE solver integrates from 0.0 to N-L The algorithm tested concluded to be quite robust. No singularities would occur. There are still, a few problems that have to be addressed which still affect the solution. Chapter 3. Creation of a Normal Directrix from Two Space Curves 3.3.2 53 Results and Present Problems Combining the equations above and integrating using the ODE class to generate the gen erators for the developable surface and the control vertices for the normal directrix yields the developable conical-like shape shown as an example in Figure 3.12. The directrix shown in the figure has actually two directrices. Both happen to coincide exactly for this example. One directrix follows exactly as the surface is being integrated. The other directrix results from solving for the normal directrix control vertices. Figure 3.12: Conic-Like Section with Flatness Tolerance Specified at 0.001 The next figure shows the same conical-like shape without a flatness tolerance speci fied. The problem of fanning is apparent at both ends. Note how the two directrices are different. The one which looks invariant is the control vertices which have been solved for and have an additional constraint of having to be perpendicuar to the end generators. The directrix which appears to follow the surface more exactly is the one which relates the apparent tangent with the current generators. The following figures display the modified conventional approach used to approximate Chapter 3. Creation of a Normal Directrix from Two Space Curves 54 Figure 3.13: Conic-Like Section With No Flatness Tolerance Specified the normal directrix control vertices which is then used by the modern approach. The examples cited are hull sections of the UBC Series fishing vessel. In figure 3.14 is an overlayed closer comparison of the two methods with the current configurations. In Figure 3.15 both versions appear to give better results. The modified approach has the tolerance set higher , 0.01, and displays the fanning problem near the stern being reduced. In figure 3.16 is an overlayed closer comparison of the two methods with the current configurations. In Figure 3.17 both versions appear to give even better results. The modified approach has the tolerance set higher , 0.1, and displays the fanning problem near the stern is almost eliminated. In figure 3.18 an overlayed closer comparison of the two methods with the current configurations is shown. Chapter 3. Creation of a Normal Directrix from Two Space Curves (a) UBC Series Surface (main deck chine and chine 1) body plan (b) UBC Series Surface (main deck chine and chine 2) profile Figure 3.14: Developable Surface Procedure Comparison, tolerance 0.001 55 Chapter 3. Creation of a Normal Directrix from Two Space Curves (a) UBC Series Surface (main deck chine and chine 1) modified approach body plan (b) UBC Series Surface (main deck chine and chine 2) modified approach profile (c) UBC Series Surface (main deck chine and chine 1) non-optimized body pian (d) UBC Series Surface (main deck chine and chine 2) non-optimized profile Figure 3.15: Developable Surface success, tolerance 0.01 56 Chapter 3. Creation of a Normal Directrix from Two Space Curves (a) UBC Series Surface (main deck chine and chine 1) body pian (b) UBC Series Surface (main deck chine and chine 2) profile Figure 3.16: Developable Surface Procedure Comparison, tolerance 0.01 57 Chapter 3. Creation of a Normal Directrix from Two Space Curves (a) UBC Series Surface (main deck chine and chine 1) modified approach body plan (b) UBC Series Surface (main deck chine and chine 2) modified approach profile (c) UBC Series Surface (main deck chine and chine 1) non-optimized body plan (d) UBC Series Surface (main deck chine and chine 2) non-optimized profile Figure 3.17: Developable Surface success, tolerance 0.1 58 Chapter 3. Creation of a Normal Directrix from Two Space Curves (a) UBC Series Surface (main deck chine and chine 1) body pian (b) UBC Series Surface (main deck chine and chine 2) profile Figure 3.18: Developable Surface Procedure Comparison, tolerance 0.1 59 Chapter 4 Optimization of a Normal Directrix 4.1 Objective Once the modern approach was deemed fairly stable the next criterion were then ap proached and listed as follows: • To better match a developable surface defined by a normal directrix to a pair of space curves. • To create an optimization function which is the mean square distance from each space curve to the nearest point on the surface. • To derive an equation to evaluate the mean square distance. • To utilize a robust non-linear optimization technique to minimize the integrated mean square distance. 4.2 Optimization Function Calculation of the normal distance between a point on a space curve and the devel opable surface, see Figure 4.1, can be obtained by three definitive relationships. This will be described and derivation of the first form of solution will be shown below: • The first relation is that the normal, n(t), of the surface is equal to the cross product 60 Chapter 4. Optimization of a Normal Directrix 61 mace curves De Jves to surface developable surface dlrectrix Figure 4.1: Relating the Distance Between Space Curves and Surface of the tangent, dp and the generator, g(t). -- - dp n(t)=xg(t) • The second relation states that any position on the surface can be calculated, q(t, .s) by moving along the directrix, p(t) and then along the generator, g(t): q(t,s) =p(t) +sg(t) • The third relation builds on the previous statement by further stating that the closest point from the surface to a point in space would be a specified distance, 1, Chapter 4. Optimization of a Normal Directrix 62 normal to the surface along the normal vector, n(t): r(u) = q(t,s) + ln(t) Substituting the second relation into the third yields: r(u) =p(t)+sg(t)+ln(t) The third relation above is the combined total of the previous two and can now be differentiated with respect to the independent variable t in order to create a differential equation which can solved. Differentiating yields the following: drdu —ds —dl —g(t)—n(t) dp = dg dn (4.1) We also have two other known relations which can be substituted, namely: - fU __.\ — 1 d p dt di (42) dt (4.3) = From these three relations we can solve for three equations and three unknowns, with the unknowns being the following dependent differentials: • The differential dependent parametric position along a space curve du dt Chapter 4. Optimization of a Normal Directrix 63 • The differential scalar dependent position along the generator g(t) ds di • The differential scalar dependent position along the normal n(t) dl di We then solve for the following relation: dr dux —g —n dr duy —gy fly dr duz —g —n 4E cit -i-- dt. dt cit dty dty -‘-i dty dt dtz dtz dtz dt (4.4) The second form involves using the three definitive relations and continuing to find relations yielding direct solutions without solving simultaneous equations. A more detailed derivation for both of the forms can be found in Appendix H. Starting with the differential relation: drdu - ds dl g- n = dt di di Taking the dot product with vector g: ( Jdu g) — ds g • g dl — g (4.6) g+sg+lg = This simplifies to the following: ds di (4.5) (4.7) (4.8) ( du (4.9) Chapter 4. Optimization of a Normal Directrix 64 Similarly, dot product with n, (4.10) and simplifying yields: (4.11) dl dt — (dr ..idu = I—n I— du dt (4.12) Similarly, dot product withy (Edu (4.13) - dt) dt — dt dt S dt dt di di Simplifying yields: (4.15) i -. . (s g + 1 n) — du di 4.2.1 414 (du 4.16 di Integral Chosen to Minimize and Simplex Parameters Used The key parameter chosen to calculate, the distance from a point on a space curve to the closest point on a surface, has to be further related to a function which describes the entire surface. Calculation of a least squares approximation for the area under the curves was chosen. One can refer to figure 4.2 for a visual representation of the area under the curves. By collecting a history of distances between the space curves and the surface, one can calculate the total area. Ideally, this would be zero. An equation in this form would work with the Downhill Simplex Method of Multidimension. The distance along the curve, the arc length, was calculated as follows: The arclength, was calculated for both curves in calculations. = jbdr (4.17) Chapter 4. Optimization of a Normal Directrix 65 Alec O.xve Und developable surface dlrecfrlx Figure 4.2: Area calculated under space curves b = I Iv(t)Idt f v(t)dt (4.18) Ja & = dS dt (4.19) Ja d pt Ja pt d = dtfa dS = v(t)dt v(r)dr=v(t) (4.20) Iv(r)Idr = Iv(t)I (4.21) (4.22) Let h 1 and h 2 be arc lengths for the two space curves r 1 and r 2 respectively. The resulting least squares type approximation for the total area would be: Area = J 1I F2 Idh 2 [1j+l2 ] 2 d1 dt (4.23) Chapter 4. Optimization of a Normal Directrix 66 From equation 4.23 the area between the space curves and the developable sur face is the function to be minimized. This is accomplished by using the equations for , 4, and for each space curve. These totalled together sum to be the degrees of freedom for the non-linear optimization method model. The non-linear optimization method used is the Downhill Simplex method. This is discussed later in this chapter. By solving all the differential equations’ independent and dependent variables from the start of the developable surface to its end, the non-linear optimization function can attempt to minimize equation 4.23. Another key point is that the dependent and independent variables sum to give ten degrees of freedom. The nonlinear optimization method also needs ten different solutions to start with. This problem was solved by moving one control vertex degree of freedom from the directrix by unit length. From these initialization procedures the Downhill Simplex method was able to start solving these equations. 4.2.2 Present Problems Both approaches stated above describing solving for the dependent variables will pro duce solutions for surfaces which do not have major changes in the rate of change of the curvature and the direction of the curvature. If the rate of change of the out-of - plane curvature with or without contribution by the in-plane curvature is very large, and changing in direction, the modern approach may have difficulty in the alignment process when trying to conform to the surface. Cases may arise, as shown in figure 4.3, where a corresponding normal on the surface which should match a position on the space curve will not occur. This may be due to the constraint that both the space curves and the surface are of finite length and the corresponding normal on the surface will point to an out-of-bounds position. One other possibility is that a localized change in curvature orients itself midway during the iteration process and another normal of infinite length Chapter 4. Optimization of a Normal Directrix 67 occurs. The first method of matching the modern approach to two space curves can also fail at a prior stage when calculating the inverse of the matrix. Reasons for a non - solution are generally for the same reasons. I DIsonce L frornspocel eto normal Figure 4.3: Possible inherent failure due to certain geometry 4.3 Downhill Simplex Method Many non-linear optimization methods as cited by Jacoby [14] would be possible to implement, but, a more suitable and generally more stable and robust method, the Simplex method [16], was chosen. A major reason for this preference was the description cited in Press [20], which briefly described its multidimensional capablities. Chapter 4. Optimization of a Normal Directrix 4.3.1 68 Explanation of The Downhill Simplex Method A brief description of the Downhill Simplex Method in Multidimensions is presented here to help clarify its use. A simplex is the geometrical figure consisting, in N dimensions, of N + 1 points (or vertices) and all their interconilecting line segments, polygonal faces, etc. [20] In three dimensions, it is a tetrahedron, as can be seen in figure 4.4. Downhill Simplex Method (a) )C (b) (C) Figure 4.4: Analogy of a Simplex For The Downhill Simplex Method The downhill simplex method has three parameters which define the moving behaviour of the simplex. These parameters are mentioned here because changes in their values Chapter 4. Optimization of a Normal Directrix 69 may be necessary. The three parameters are defined by Nelder [16] as a, 6, and . These parameters correspond to: reflection, contraction and expansion. The reflection parameter, a, is a positive constant, which is a scalar multiplication constant which mirrors the point through the simplex. The contraction parameter, 3, is a constant which values lie between 0 and 1. It is a ratio of the distance of the point relative to the simplex centroid. The expansion coefficient, , is greater than unity and it is a ratio of the current point to the centroid with a point along the line joining the point to the centroid. For a more detailed explanation refer to Nelder [16]. For detailed coded form of the Downhill Simplex Method, programming and tests performed, refer to Appendix L. 4.3.2 Results and Present Problems As stated previously, surfaces which resemble close to developable surfaces to begin with will be more likely to have a solution that closely matches the original space curves. Below are two examples which demonstrate both simple and difficult solutions encountered. 4.4 Examples A few examples are shown in this section to help explain the various phenomenon of the results of the theory created this far in the research. 4.4.1 Example of Single Surface with Conical Properties The first example demonstrates how the non-linear optimization can successfully solve a free-form surface with conical properties. Shown below in figure 4.5 are two views of a conical like surface. Note that due to the symmetry of the object in figure 4.5 the generators also follow symmetrically along the surface. Chapter 4. Optimization of a Normal Directrix (a) Conical-like developable surface view 1 70 (b) Conical-like developable surface view 2 Figure 4.5: Conical-like developable surface 4.4.2 Example of UBC Series Demonstrating Present Problems The UBC Series vessel was posed as the immediate goal to make into a developable ship series. In figure 4.6 the surface tested was using the 1st and 2nd chine as the space curves. As seen in figure 4.6 the optimized version had difficulty with the adjustments due to the in-and-out-of-plane curvature that results from this type of surface. Referring to figure 4.7 both versions appear to give moderately successful results. The modified approach has the tolerance set low near the stern. , 0.001, and displays the fanning problem Chapter 4. Optimization of a Normal Directrix (a) UBC Series Surface (Chine 1 and 2) nonoptimized (b) UBC Series Surface (Chine 1 and 2) opti mized (c) UBC Series Surface (Chine 1 and 2) com parison Figure 4.6: Developable Surface Failure 71 Chapter 4. Optimization of a Normal Directrix (a) UBC Series Surface (main and chine 1) modified approach body plan (b) UBC Series Surface (main and chine 2) modified approach profile (c) UBC Series Surface (main and chine 1) optimized body plan (d) UBC Series Surface (main and chine 2) optimized profile Figure 4.7: Developable Surface success, tolerance 0.001 72 Chapter 5 Intersection of Developable Surfaces and Flat Plate Layout Once developable surfaces have been created, it is necessary to define the edges of those surfaces by intersection with adjacent surfaces. An additional consideration is once the bounded developable surface has been created a flat plate layout is necessary. A derivation of the solution is also presented here. 5.1 Intersection of Developable Surfaces Once the developable surfaces have been created there is the requirement that they must intersect without gaps or spaces between them. Shown below is the vector representation for initial conditions of the surfaces and their orientation relative to each other. In figure 5.1 we see the various vectors that are dependent upon each other. We can show the vectors dependence on each other as follows: Figure 5.1 presents the independent parametric variable I. The two dependent para metric variables are u and v. The independent space curve is p(t); one dependent space curve which is a function of u is f(u); the other dependent space curve is g(v), which is a function of v. The curves p(t), f(u) and g(v) are normal directrices representing the developable surface of interest, and the adjacent surfaces on each side. The intersection curves joining two developable surfaces are r 1 which is relating f(u) and p(t) and the other intersection curve r 2 relates g(v) and p(t). A detailed derivation of the equation used to calculate the intersection of the devel opable surfaces are located in Appendix I 73 Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 74 A 4 gft: 72 p(\X*/ S b(k ...• dh h(v) Figure 5.1: Intersection of Three Developable Surfaces 5.1.1 Derived Equations Yielding Intersections of Surfaces Derivation of the equations used to calculate the intersection of the developable surfaces are as follows: 1 r = Pt + = dp dg 1 ds +sit--+--gt 1 dr —k- it (5.1) (5.2) Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout (E + •. 75 (5.3) dt gt— - dg dt gt’dp - 4. - du Igu S df — (5.4) dp d) — S Idfu (5.5) df (5.6) drdu (5.7) - 1 dr -;u- = But, s2ugu and, r 1 duQfu u + S2 dg + fJ du — (5.9) = = = f + s2ugu, and — (5.10) f (5.11) p+sg (5.12) Substituting gives, dfdf (5.8) (5.13) = f 2 d _S2u) dt u du du(dhdh (5.14) Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout / —s 76 __.. (dp +it+- du du_ —— dt (5.15) ‘— 1df—.——-..(pt+sitgt df f 2 d ---i -- dt = 2 r — Pt + --- / —k dp dp (5.17) S2tgt - dgt + + _ .... dgt • 2 r (dIi Ija x g) I____. \ YtIdp dg dt /____. (dh (5.19) /___. g (5.18) +S2r !__• cit gt 2 ds —— — (d1i ã._xv) — (5.16) du dt — - 2 dr _fu)) (5.20) p 2 Id = E\ h+ (5.21) (5.22) Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout dr 2 dt dv 2 dr — 523 dvdt — 2 dr (5.24) = 2 r 2 r 2 dh dr (dh h+s 1 (5.25) 2 r (5.26) — p+ h 2 d tdt 2 J d1z h 2 d dgdgdv dtdvdt 5.1.2 (5.27) sä dv (dh dh = ( Idt dv 77 (5.28) dv —1i) 530 Results and Present Problems The intersection of developable surfaces was tested and proved to be fairly stable when well behaved surfaces were created. Conical Type Surface Tested An example shown in figure 5.2 displays very tightly intersecting surface without any problems. The intersection algorithm needs N + 2 surfaces when designing for N surfaces. In the next subsection the algorithm is shown using what we term phantom surfaces. Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 78 IC I \\\\\\\ (a) Conical-type surfaces view point (1,1,1) (b) Conical-type surfaces front view Figure 5.2: Conical-type Surfaces Intersections Criterion of Phantom Surfaces The criterion of phantom surfaces was introduced to provide a secure mechanism for preventing interconnecting surfaces from having “cracks” or “leaks”. The conical-type surfaces presented in figure 5.2 show surfaces which have been created using the modified approach. These all appear to be developable but only the inner two have been con firmed as developable using the modern approach. The two outside surfaces were used as phantom surfaces in order to provide an intersecting edge that the inner two could use. Looking at the figure 5.3c) and figure 5.3d), these two surfaces use the modern approach and yield aligned surfaces with the intersection algorithms. The cracks will be more evident when looking at the ARCTIC series of fishing vessel or the UBC series vessel. Figure 5.3c) and figure 5.3d) are shown to display the phantom surfaces Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 80 Flat Plate Layout 5.2 In order to convey the information of a developable surface into a useful form, the fiatplate-layout was created so the user can essentially use this output as a template. 5.2.1 Derived Equations giving Flat Plates Dimensions and Interplate An gles Once a developable surface has been created, the shape and position of the consecutive segments that make up the developable surface need to be known. Below is the derivation of the algorithm used to provide the information of the shape of the segments in the x y plane and the angle and rate of change of bending of the segments relative to each previous one. For a more detailed explanation of the derivation one should refer to Appendix J / !2 dt dt Lht Figure 5.4: Flat plate layout derivation The resulting differential equation is used to find the new tangent and corresponding Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 81 point on the x-y plane where the next point of the plate is located: a12 q \dt2 dt) dt UP (5.31) = 5.2.2 Format of Output The example of the flat plate layout for the developable mobius strip can be seen in figure 5.5 (a) 3-D view of developable mobius strip (b) Corresponding flat plate layout Figure 5.5: Developable surface and flat plate layout The flat plate layout corresponding to the numbered plates in figure 5.5 has data listed in table 5.1. Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout Plate no. Parametric value out of plane angle (deg) 0 0.00 4.641389 0.975138 1 0.25 13.924170 8.451800 2 0.50 23.206951 19.508427 3 0.75 32.489716 27.889669 4 1.00 32.115269 21.843246 5 1.25 29.928911 13.160947 6 1.50 38.130260 19.531147 7 1.75 56.255440 60.359379 8 2.00 55.513470 41.122555 9 2.25 34.989170 8.168997 10 2.50 21.341158 3.404501 11 2.75 14.529543 2.682558 12 3.00 12.068573 6.007172 13 3.25 9.672798 13.691539 14 3.50 10.510071 9.681314 15 3.75 15.940498 9.440858 16 4.00 18.051476 19.695953 17 4.25 12.895594 21.226387 18 4.50 7.738149 11.441552 19 4.75 2.579468 4.393825 20 5.00 -0.016314 0.398113 Table 5.1: Flat plate data 82 Chapter 6 Choice of Computer Language and CAD Program When this thesis work was initiated, one of the objectives was to initially create generic tools in a very portable modular language and then apply them to a well established open ended CAD package that is used throughout industry. At the present stage of the work both objectives seem to have been met. 6.1 Computer Language Chosen The choice of which computer language to program was based on a number of design constraints. The language chosen should be flexible in case the method of design had to be altered. The three designs considered were: • Top-down structured design • Data-driven design • Object-oriented design Of these three types of design, data-driven design was initially dropped because of the nature of the design tool. Because of the financial constraints on the project and the software languages available on the platform chosen, the top-down structured design was initiated. This was then taylored into a modular top-down structured design. After the release of a popularly supported compiler for the platform was chosen, the design was then switched to an object-oriented design. The language initially chosen was ANSI C, best supported at the initial time of development by Microsoft for the PC platform. 83 Chapter 6. Choice of Computer Language and CAD Program 84 Later, after an initial attempt with Walter Bright’s, Zortech C++, the switch was made to Borland’s C++, which proved to have excellent support tools, a good development environment and better technical support. 6.1.1 OOP - Object Oriented Programming When it was realized that a true portable object-oriented language was available, the design was taylored using influences from OOP, object oriented programming, OOD, object oriented design, and OOA, object oriented analysis. First a definition of each is in order: • OOP Object-Oriented Programming Object-oriented programming is a method of implementation in which programs are organized as cooperative collections of objects, each of which represents an instance of some class, and whose classes are all members of a hierarchy of classes united via inheritance relationships [5]. • OOD Object-Oriented Design Object-oriented design is a method of design encompassing the process of object oriented decomposition and a notation for depicting both logical and physical as well as static and dynamic models of the system under design [5]. • OOA Object-Oriented Analysis Object-oriented analysis is a method of analysis that examines requirements from the perspective of the classes and objects found in the vocabulary of the problem domain [5) Chapter 6. Choice of Computer Language and CAD Program 6.1.2 85 Portability to Different Platforms One major consideration when doing software development for an application is the standardization and portability of the code being developed. Up until a few years ago, one of the most standardized and portable languages was ANSI C. Now with more features and better type checking, ANSI C++ appears to be the successor. Porting some of the sample computer code mentioned below to different compilers and to different computer platforms proved to be a very simple and straight forward task. Very little problems were encountered. The sample compilers tested were Borland C++ 3.1 and Zortech C++ 3.0. The platforms tested were DOS based PC platforms and the SUN SPARC UNIX platform. No major problems were encountered. 6.1.3 C++: Classes and Features The language chosen to best implement the work being done was the C++ language. This language proved to be very much more than just a superset to the C language. Many of the class features, when used in its entirety, better shape the programming environment into a true object oriented program and design. A little explanation is presented in the next sections as to how the tools were used to better inform the reader as to the approach taken. Class vector *iaclude <iostream.h> *ifndef _VECTOR_ *define _VECTOR_ en direction class matrix; { I, Y, Z }; Chapter 6. Choice of Computer Language and CAD Program class vector mt 86 { n; static jut default_length; float *eleent; public: vector 0 {ndefault _length; ele.entnev float En ; }; vector(int din) {ndi;eleentnew float[n];}; vector(int di., float x); vector( coast vectorl v vector0; static void set_default(jut n){default_length=n;}; static void set_default 0{default_length3;}; float& operator[J( mt 1) {return ele.ent[i];}; float operator[J(int i) coast {return eleent[i] ;}; friend mt size(const vectort v) { return v.u;}; vectort operator=(float x); vector& operator(const vectork v); vector& operator(const .atrix& a); friend vector operator+(const vectori a, coast vector& b); vector& operator+( coast vectork a); friend vector operator—(const vectork a, const vector& b); vectorl operator—(const vector& a); friend vector operator*(const vector& v, float a); friend vector operator*(float a, coast vector v); friend float operator*(const vector& a, const vectort b); II dot product operator friend vector operator*(const matrixk a, coast vectork b); vector& operator*(float a); vectort operator*(const vector& a); II e1e.entby—element nltiplication friend vector operator/(const vectork v, float a); friend istrea.& operator>>(istreamk, vector&); friend ostrea.k operator<<(ostreamk, const vector&); *endif The class vector is the name of a class which enables the user to use N dimensional vectors for calculations. A C++ feature, termed operator overloading enables the user to write equations almost exactly as one would write them by hand. Chapter 6. Choice of Computer Language and CAD Program 87 Included below is the vector.h class header file which shows how the class is initially designed. Some explanation is in order at this stage to clarify to the reader these basic features which are used throughout the development of other classes used in developing the various design tools. Looking at the listing of the class vector.h, the compiler directives, #ifndef _VECTOR_ and #define NECTOR_ are compiler directives which will not enable redeclarations. Be low is an optional enum direction { X, Y, Z }; which enables the user to declare X as a 0, Yas land Zas 2. The next line class matrix; declares the matrix class before class vector.h can be declared, thereby, enabling class matrix to be used inside the vector.h class. On the next line below is the class declaration class vector {, which assigns a class name. The next line, since it does not have any member access control stated as private:, protected:, or public:, then the members of a class declared with the keyword class are private: by default. The next line, whose member access control is public: shows some very interesting features of C++. The first constructor, depicted by vector() (n=3; element=new dou ble[n]);, which shows a default size, n, for the class to allocate memory. vector() is an operation which is empty of a type. This case is where the type is hidden and some mechanism must be provided for a user to initialize variables of that type. Located on the next line is, vector(int dim) n=dim; element=new double[n], which is a constructor of type double. Moving down one line, vector(const vector & v), is a copy constructor. “A copy constructor for a class vector is a constructor that can be called to copy an object of class vector that is, one that can be called with a single argument of type vector. A copy constructor is generated only if no copy constructor is declared”. The next line, iector() delete element;;, is contains what termed a destructor. A destructor automati cally deallocates dynamic memory when a class goes out of scope. Moving to the next Chapter 6. Choice of Computer Language and CAD Program 88 line contains one possible case of operator overloading. The line, doubleF_4 operator[](int i)returri element[i];;, enables the programmer to use the [] operator with the vector class. For example, the following code fragment would look like: main() { vector x; x[0] = 0.0; x[1] = 1.0; x[2] = 2.0; cout<<x[0] <<“ x[1] << “<< “ “ << x[2] <<‘\n’; } There is also another form of the operator [] which implements the named const. There is also the following code fragment which resides inside of class vector. It is in the following form, float operator[] (mt i) const{returri element[i];};. This is necessary since the compiler complains when several operations are done successively and invoke temporary storage variables. For example: vector x,y,z,a,b,c,d; .x,y,z,a,b,c,d . . .get assigned values Chapter 6. Choice of Computer Language and CAD Program 89 then x (y*z) *a* (b*c) *d; cout<<x; The above list of code shows how the operator jj is used as well as the << stream operator is overloaded to standard output for the class vector. Syntactically this code is much simpler and far more legible to the reader. This code is simple and easy to follow thereby leaving the details of implementation to the compiler. Another note one should make is that the declarations of mt n; and double *element; are private inside class vector. Therfore, any queries of these variables must be done by a member function within classvector. So, for example, in order to find the size declared class of type vector, one has to invoke either friend v.n;; or mt ii of a size (vector éY v)return mt size ()return n;; in order to access items from the private section of the class. The next three lines of code depict three types of operator overloading one could encounter when developing classes. The first operator, vectorJ operator_—(vectoré4 v);, equates two classes of type vector, this member function also utilizes the *this pointer, which is a special reserved pointer used in C++. The details of how the information is passed is left for the reader to investigate. The second operator, vectors operator=(const matrix& a);, equates two objects of class type vector and matrix. This enables an equat ing of classes of different types. The next example operator,friend vector operator *(const vectors a, const vectors b);, requires a friend to be declared. This code is optional since the same feature could be invoked by the following, vectors operator*(vector a);. Chapter 6. Choice of Computer Language and CAD Program Class matrix In the listing below, the matriz.h class header file implements the vector.h class. tinclude <vector .h> tinclude <iostream.h> #ifndef ..AATRIX_ idefine _JIATRII_ class matrix { tnt nr, nc; float **element; public: matrix(int rows, mt columns ); matrix(const matrixi a); matrix(const vectork a); matrixfl; matrix& operator(float x); matrixi operator(const matrix*); float* operatorfl (tnt i) {return element[i] ;}; coast floata operatorS (hit i) coast {return element[i] ;}; vector row(int i); friend 1st n..rows(const matrixi a){return a.nr;}; friend tnt n_coluns(const matrixt a){return a.nc;}; vector column(int i); friend matrix exp(const matrixt a); friend matrix element_mult(coast matrixi a, const matrixi b); friend vector diagonal(const matrixi a); friend matrix operator+(const matrix& a, coast matrix tb); friend matrix operator-(const matrixt a, coast matrixi b); friend matrix operatore(const matrixi a, coast matrixi b); friend vector operator*(const matrixt a, coast vectort b); friend matrix operators(const matrixt a, float b); friend matrix operatox*(float a, coast matxix& b) {retunt baa;); friend matrix transpose(coast matrixk a); friend matrix solve(const matrixi a, coast matrixt b, float adet friend matrix identity(int); friend matrix inverse(coast matrix& a, float* detNULL); friend istreamk operator>>(istreama, matrixi); N1JLL); 90 Chapter 6. Choice of Computer Language and CAD Program 91 friend ostreant operator<<(ostreaak, const satrixa); friend hit operator<(const .atrix&, conat •atrixt); friend zatrix abs(const aatrix&); I; *endif One can see the similarity between class vector and class matrix and how they can interact. Class Curve and Surface tinclude “matrix .h” tinclude <fstream .hpp> *ifndef _CURVE_ tdefine _CURVE_ class Curve{ mt n; char splinetype; vector epoint; static matrix SB; public: void sptype (char what) {splinetypewhat ; char typeofsplineO{return splinetype;}; int sizeO{return n;}; CurveO{n5; splinetype ‘1’; point = new vector[51;}; Curve(int nua){nnum; splinetype’l’; point Curve(Curvek a); CurveO{delete[n] point;}; void realoc(int na); vectork operator El (hit i){return point Curvet operator(Curvet c); vector operatorO(double t); vector tangent(double t); vector curvature(double t); double deriv_2c(double t, hit i); double deriv_c(double t, mt i); Ii) ;}; = new vector[n];}; Chapter 6. Choice of Computer Language and CAD Program static void initbasis(char sptype’2’, double b11.0, double b20.0); static atrix BEB; tendif tifndef _SURt *define _SURF_ class Surface{ private: Curve *leftedge, erightedge; Curve *directrix; Surface *leftsurf, *rightsurf; public: Surface( Curve *leftedge, Curve trightedge, Surface *leftsurf, Surface *rightsurf,int ncontrolpntss,double step7.0, double outpO.0S,iut surfnol); void Rightsiuit(Surface trights); void Leftsinit(Surface Clefts); double Optiaize(double toleranceo.2, mt itermax = 250); void Develop(void); void Developt (void); void Developtl(void); void Drav_3d(int gemtoso,double step7.0,int surfnol); void Draw_33d(int gennosO,double stepfl.O,int surfnol); void Draw..3id(int genno5O,double stepfl.0,int surfnol); void Drav.2d(int gennoso,double steopfl.0,int surfnol); SurfaceO; Surface(Surfacet s); SurfaceR operator(Surface& s); 1; tendif *ifndef _FILELISTS_ Idefine _FILELISTS_ tinclude <fstrea. hpp> class Filelisto{ public: 92 Chapter 6. Choice of Computer Language and CAD Program 93 ofstrea. file; Filelisto *next; ofstrea.k operator[](iut n); class Filelisti{ public: ifstream file; Filelisti *next; ifstream& operatorlj (jut a); *endif The following classes, Curve arid Surface also implement vector and Curve. These sort of implementations were able to make the code more legible and apply the testing of Developable Surface modules much earlier in the design cycle. Class ODE *ifndef _ODE_ *define _ODE_ *iaclude “vector .h” class dynamic_system { private: double time, step_size; vector state; vector error_scale; vector (ederivative) (double, vector&); public: dynamic_system(vector (*)(double, vector&), double start_time, vectork initial_state, vector *error_scaleO); double& whenO; void resetO; double& operatorO(int i); vector operator() (double t); vector step(double delta); Chapter 6. Choice of Computer Language and CAD Program 94 ‘- = I I I I I \\ Figure 6.1: Computer Platform Selected Initially void rk(double tO, double& del_t, double ti, vectorl x, vector (sf)(double t, vectort x), vectork err_scale); #endif The above code class ode, ordinary differential equation solver, also shows another useful combination of class building to contribute to a library of useful class tools. 6.2 Computer Platform Selected Chapter 6. Choice of Computer Language and CAD Program 6.2.1 95 Pc 386/486 Environment The development of the theory for this thesis and the validation via progranmiing was implemented on a PC. The PC, personal computer, platform was decided upon because it is the most popular and inexpensive platform in use for almost every application today. The PC is primarily an opened ended architecture in which a wide variety of hardware can be interfaced. With PCs becoming faster and less expensive every year most of the speed constraints of the thesis programs were relaxed. At the present state of the thesis researçth and design, the programs execute fast enough to become interactive. Commercial viability should be considered. 6.2.2 computer Language Selection in this Platform At the start of this thesis the computer language chosen was the C programming language. This language proved to be sufficient for initial implementations. Further development of computer languages, such as C++, proved to offer much more sophisticated features and modularity of code that would offer more efficient development cycles. This proved correct and the theory, code, test, development cycle time was greatly reduced. Many other features of the C+++ language also proved invaluable towards imple menting the theory. One of the most used features of the C++ language was the use of operator overloading. This enabled vector calculus equations to be written in code in almost the same convention one would use when writing by hand. Other features which were discussed above all contributed to more legible code and as a collection of useful tools which other applications would ensue. Chapter 6. Choice of Computer Language and CAD Program 96 CAD Program Selected 6.3 During the course of development of the theory a decision was made as to whether or not writing a computer platform independent CAD package should be undertaken. After assessing the existing CAD packages used on several computer platforms, the conclusion was that one would be wiser to develop the special code to an existing CAD package. The CAD package chosen was assessed as the most popular and efficient to develop upon. This CAD package was AutoCAD. During the course of development of the theory, ongoing assessment was done on AutoCAD’s market position and development tools. This proved to be benefitial to this project as well as it reflected upon us how strong a commercial product this was in the CAD market as well as their desire to move to other hardware platforms. AutoCAD was reviewed and concluded to be the most popular CAD package in in dustry and contained the best development design tools out of all of the CAD packages surveyed. Initial development started with AutoCAD release 10 on the 386 personal com puter. This proved to be frustrating and fell short of our expectations and requirements. The release of AutoCAD release 11 and 12 proved to answer most of our questions and solve most of our earlier problems. Now pending publishing of this thesis ongoing devel opment is being done towards commercialization of developable surfaces as a third party developer package. Much effort is needed in order to create user interactive tools sophis ticated enough to ensure useful interactive CAD tools that the user would find intuitive to use. 6.3.1 CAD Program Environment and Open Architecture The CAD program environment using AutoCAD ADS development tools proved to be what was needed for developable surfaces. The ADS development tools are C language Chapter 6. Choice of Computer Language and CAD Program 97 calls which are registered and loaded as applications for AutoCAD. The ADS develop ment tools have proved to be of high quality and as an open architecture to build and expand upon the existing tools contained in AutoCAD. This open architecture has en abled AutoCAD features and newer releases to grow very rapidly. The writing of ADS development tools in ANSI C language has also enabled AutoCAD to be ported to other platforms very quickly and easily with minimal effort when compared to other languages. 6.3.2 Present Limitations So far the present limitations of AutoCAD and ADS development tools are few. Some awk wardness has occured when implemented AutoCAD ADS applications tools with some of the C++ modules. True C++ OOD, OOP techniques can be affected by some of the C functions and how the ADS environment is structured which limits some of the C++ design. Future releases of AutoCAD and ADS tools will probably evolve into C++ classes and methods. Existing design changes are being implemented quite easily into the Auto CAD ADS environment. Autodesk has been very supportive and provided encouragement for our ongoing research and development implementing C++ to ADS. Chapter 7 Demonstration Examples In this chapter there are four example sections which display the benefits and problems still facing the present design. The first section describes the power and flexibility of the base module which, when unconstrained, can create more difficult geometric anomalies in computational geometry. The other sections show more realistic applications and the difficulties that arise trying to “best fit” two or more space curves to the exact geometry desired. One must be made aware that from any two or more space curves, a developable surface may not exist; only a closest fitting one may be a solution. The reader should be made aware that this is a design tool, and it should also provide the user with feedback as to what type of geometries are potentially developable. 7.1 Developable Mobius Strip The first example figures shown below in figure 7.1, which were presented in chapter 3, show the flexibility of the base module, given freedom to solve the geometric problem when working from a directrix and two generators only. This geometric anomaly is interesting to view and cite as an example and perhaps as mathematical art, but provides no potential commercial benefit to the industrial designer. The next few sections exemplify more practical examples in commercial areas such as manufacturing. When “harder” constraints are involved new problems arise. Examples could have been cited where these design algorithms would prove to the reader that they 98 Chapter 7. Demonstration Examples (a) Developable mobius strip view 1 99 (b) Developable mobius strip view 2 Figure 7.1: Developable mobius Strip Chapter 7. Demonstration Examples 100 work very well, but, they would also prove to be misleading. This thesis was presented in a very objective engineering manner which shows both the positive and negative sides inherent in design. 7.2 Simple Conical Developable A simple conical developable surface was created in order to verify and check compliance with known developable surfaces, fanning and intersection of surfaces. Shown above 7.2 are surfaces which intersect and verify results. 7.3 Arctic Fishing Vessel 7.4 UBC Series Fishing Vessel Referring to figure 7.5 the ubc series was found to be the most difficult to model since the bow contains compound curvature hull form in certain regions from station 0 to station 2.0. Chapter 7. Demonstration Examples 101 (a) Conical-type surfaces view point (1,1,1) (b) Conical-type surfaces front view (c) Conical-type surfaces view point (1,1,1) (d) Conical-type surfaces front view Figure 7.2: Conical-type Surfaces Intersections Chapter 7. Demonstration Examples 102 (a) ARCTIC vessel view point (1,0,0) (b) ARCTIC vessel view point (0,4,0) (c) ARCTIC vessel view point (0,0,1) (d) ARCTIC vessel view point (-1,-1,0.2) Figure 7.3: Arctic Vessel Conventional Approach Chapter 7. Demonstration Examples (a) ARCTIC vessel view point (1,0,0) J// 103 (b) ARCTIC vessel view point (0,-1,0) ‘I’ll! \\\\‘\‘\‘ \ (c) ARCTIC vessel view point (0,0,1) (d) ARCTIC vessel view point (-1,-1,0.2) Figure 7.4: Arctic Vessel Modern Approach Chapter 7. Demonstration Examples 104 (a) UBC series vessel view point (1,0,0) (b) UBC series vessel view point (0,-1,0) (c) UBC series vessel view point (0,0,1) (d) UBC series vessel view point (-1,-1,0.2) Figure 7.5: UBC series Vessel Conventional Approach Chapter 8 Conclusions and Recommendations 8.1 Conclusions The research objectives were • to develop an algorithm in order to find a normal directrix such that the resulting developable surface lay close to two space curves representing desired edges of the developable surface. • to create an algorithm to intersect developable surfaces and to generate the flat plate layouts and angles. • to implement these algorithms using a modern computer language and a popular CAD package in order to assess the practicality of the approach. Normal Directrix from Two Space Curves An algorithm was created to compute a normal directrix for a developable surface from a pair of space curves. Differential equations were derived for defining a set of generators lying between the space curves, similar to the approaches taken by Nolan and Clements. The set of generators were used to compute a normal directrix, which was approximated by a parametric spline. To ensure that the rulings between the ends of the space curves were also generators of the developable surface, the tangents at the ends of the space curves were aligned so that the cross-products between the end rulings and the two tangents at each end 105 Chapter 8. Conclusions and Recommendations 106 were co-linear. The control verticies of the normal directrix were adjusted so that the developable surface started from one end of the normal directrix aligned with the opposite end generator. The algorithm created yielded undesireable fluctuations in generator direction when the surface was nearly flat. A threshold limit for the curvature was used, below which the generator direction was held constant. Practical examples with hard-chine ship hulls demonstrated the robustness of the approach. A non-linear optimization technique, the downhill simplex method, was implemented to further refine the shape of the developable surface. Limited success was achieved. Intersection of Developable Surfaces and Flat Plate Layout An algorithm was derived to define the edges of a developable surface by intersection with adjacent developable surfaces. The method proved to be robust, with the exception when two adjacent developable surfaces are nearly co-planar. Differential equations were also derived for creating a flat plate layout of a developable surface, including axes of bending and plate curvatures. Implementation The C++ programming language was used to implement the algorithms. This enabled an object-oriented and modular approach, for example 3-D space curves were implemented as a class of objects with member functions such as position, tangent and curvature as a function of the independent parameter. As a class of objects, details of the implemen tation of curves were transparent to or isolated from the rest of the program. Other features of this programming language enabled the main program segments to be written clearly and concisely in a mathematical format. Chapter 8. Conclusions and Recommendations 107 AutoCAD was selected as the host CAD package because of its wide acceptance, multiple platforms and development tools. Only preliminary steps have been taken to date to incorporate the algorithms and code developed directly into the CAD package. 8.2 Recommendations • Alternate approaches to the threshold limit on curvature should be considered. • Drop non-linear optimization for adjustment of the shape of the normal directrix. • Implement a threshold limit for the shape of the intersection between adjacent surfaces when nearly co-planar. • Implement non-uniform rational beta spline and non-uniform rational tension Catmull Rom spline when they are available. • A basic set of algorithms should be implemented including — — — — creating a normal directrix from two space curves, intersection of adjacent developable surfaces, generation of flat plate layout and direct interactive control of the shape of a normal directrix. Bibliography [1] R. A. Adams. Calculus of Several Variables. Addison-Wesley Publishers, 1987. [2] Glenn D. Aguilar. Definition of Developable Surfaces with High Level Computer Graphics. In Proceedings at the Pacific Northwest Section of the Society of Naval Architects and Marine Engineers, pages 1 [3] B.A. Barsky. — 21, January 1987. Computer Graphics and Geometric Modeling Using Beta-Splines. Springer-Verlag, 1988. [41 P. E. Bezier. Numerical Control-Mathematics and Applications. John Wiley & Sons, 1972. [5] G. Booch. Object Oriented Design with Applications. Benjamin Cummings, 1991. [6] W.E. Boyce and R.C. DiPrima. Elementary Differential Equations and Boundary Value Problems. John Wiley & Sons, 1977. [7] J.C. Clemens. A Computer System to Derive Developable Hull Surfaces and Tables of Offsets. Marine Technology, 18(3):227 [8] T.D. DeRose and B.A. Barsky. — 233, July 1981. Geometric Continuity and Shape Parameters for Catmull-Rom Splines (Extended Abstract. Proceedings of Graphics Interface, (27):57 — 64, May - June 1984. [9] T.D. DeRose and B.A. Barsky. Geometric Continuity, Shape Parameters, and Ge ometric Constructions for Catmull-Rom Splines. ACM Transactions on Graphics, 7(1):1 — 41, January 1988. 108 Bibliography 109 [10] A. B. Dunwoody. Computer Aided Design of Developable Surfaces. not yet pub lished, 10, 1989. [11] editor E.V. Lewis. Principles of Naval Architecture, 2nd ed. Volume 1,2,3, Society of Naval Architects and Marine Engineers, 1989. [12] J.D. Foley, A. VanDam, Steven K. Feiner, and John F. Hughes. Computer Graphics Principles and Practice. Addison-Wesley Publishing Company, 1990. [13] M.E. Hohmeyer and B.A. Barsky. Rational Continuity: Parametric, Geometric, and Frenet Frame Continuity of Rational Curves. ACM Transactions on Graphics, 8(4):335 — 359, October 1989. [14] S.L.S. Jacoby. Iterative Methods for Nonlinear Optimization Problems. Prentice Hall, 1972. [15] U. Kilgore. Developable Hull Surfaces. Fishing News (Books) Ltd., 1967. [16] J.A. Nelder and R. Mead. A Simplex Method for Function Minimization. Computer Journal, 7:308 — 313, 1965. [17] T.J. Nolan. Computer-Aided Design of Developable Hull Surfaces. Marine Tech nology, 233 — 242, April 1971. [18] J.C. Beatty R.H. Bartels and B.A. Barsky. An Introduction to Splines for Use in Computer Graphics and Geometric Design. Morgan Kaufmann Publishers, 1987. [19] D. F. Rogers. B-Spline Curves and Surfaces for Ship Hull Definition. Society of Naval Architects and Marine Engineers, 1977. [20] S.A. Teukolsky W.H. Press, B.P. Flannery and W.T. Vetterling. Numerical Recipes in C. Cambridge University Press, 1988. Appendix A Mathematical Notation for Partial Differentiation Notation used throughout is that which is defined in the text by Adams[1. Notation example: p = f(x,y,z) where x = x(t) y = y(t) z = z(t) We can then say: p = F(t) = f(x(t),y(t),z(t)) where we can regard p as a function of the single variable t. Since p depends on t through both of the variables of f, the chain rule for three terms: dp dt — 7 a dx op dy dz + + Oxdt 9ydt Ozdt The use of the straight “ d” denotes derivatives of only one variable. For example: The Derivative • while op means 1 f ( x,y,z) rIO means F’(t) = a —f(x,y,z) has Appendix B Derivation of Developable Surface B.1 Constraints which Define the Developable Surface The definitions presented below were from previous work and on going research. .i!LtT g+.JT g g dt dt 2 dt Figure B.1: Derivation of Developable Surface Vectors Must Be Perpendicular (g+ d d T)(+LT) 0 d dg dp (+g dg dt 111 dp dt = 0) = 0) — — cPp 2 dt g Appendix B. Derivation of Developable Surface 112 Vector g is of Constant Length dg The Normal is Invariant Along a Generator b a x(r,t) g(t) r p(t) t Figure B.2: Derivation showing normal is invariant along a generator x(r, t) = p(t) + rg(t) where r is a scalar dx dp dg = b= dx -;--=g ar n = dp dg normal=x=(+r)xg = dp dg (-jjxg)+(r--xg) Appendix B. Derivation of Developable Surface But dp (- x g) therefore B.2 dg ( x g) ze. x g). are 113 parallel 0 = Proof Using Constraints dg_ dt (.g)dp (2)dt Let B be the constant we need to satisfy: dg dt —B dt CONSTRAINTS 1. Vectors must be of unit lentgth. g•=O==g•B=0 2. Normal is invariant along a generator. (gx). 0 =(gx).B = 0 Identity: A. (B x C) = B. (C x A) B5.(gx) = 0 Bg.(x) = 0 Appendix B. Derivation of Developable Surface 114 3. Vectors must be Perpendicular. dg dt B “dt dp dt — — dt’ p 2 d — — — d p 2 g 2 dt (2. _dt2 g - \dt therefore: — dt = dt g) dp (.)dt Appendix C Derivation of Rate of Rotation of Generator Differential Equation This derivation also includes some new terminology. A differential equation is derived which calculates the rate of rotation of a generator with respect to changes in one of the parameters of the directrix curve. NN’ Theta Phi Directrix Figure C.1: Vector locations and corresponding angles The angle Theta, herein referred to as 0, is the out of plane rotation of G’ with respect toG. The angle Phi, is herein referred to as , and is the in plane rotation of G’ with respect to G. The variable a is a parameter describing directrix. The rate of rotation of a generator with respect to changes in one of the parameters 115 Appendix C. Derivation of Rate of Rotation of Generator Differential Equation 116 of the directrix curve is written as: d 0 2 dadt The desired expression is the rotatior of a generator with respect to changes in one of the parameters of the directrix, which is written as: dO da This desired expression can be defined as: da where, N is da unit normal defined as: the N Txg therefore, £4 dtda --(.N dt\da — giving, d 8 2 dadt — — d dg g 2 dadt dt dN dl From the derivation of a developable surface: 2 \dt dl Vdt dt where T is of unit length: T = dt /.4R Vdt dt therefore, dadt Vdt dt da da / Vdt dt Appendix C. Derivation of Rate of Rotation of Generator Differential Equation 117 The second term in the above equation disappears because the derivative of a constant is zero. therefore, (g)dT — dadt Vdt dadt dt dt \da dt J (C.2) One of the next terms to define is the first term of the second part of = da çT + ON where, = dT g aa _—— combining gives, — (. g) T + ON (C.3) One of the next terms to define is the second term of the second part of Using the differentiation product rule on one of the definitions: N=Txg yields, dN —- /dT = ‘\ x g) + / x dg’\ (C.4) Appendix C. Derivation of Rate of Rotation of Generator Differential Equation 118 Combining Equation C.3 and Equation C4 to form the second part of equation C.1, da dt (—(.g)T+oN). ((xg)+(Tx — (4 x g)) + (_ g) (T. = . ((_z . )) g) (T. (T x (C.5) + 9(N.(xg))+6(N.(Tx)) several terms drop out of Equation C.5: 0 (N. (Txff)) = 0 since, BT, hence, T x T = 0 Using Identity A. (B x C) = B.(CxA) (4JZ(g xN)),butT = gxN dt = Another term also drops out: o(N.(x)) resulting in -- T since, T remains 0, unit length One other term drops out: (—.g) (T.(Tx)) but dt Therefore dg dN da dt Using a vector triple products dot = BT, resulting in (T x T) /dT ‘\( IdT ——g) T.—jj-xg and cross identity = 0 Appendix C. Derivation of Rate of Rotation of Generator Differential Equation IdT ‘\ IdT -_.g)-.(gXT) = using the fact that N Txg fdT da = = ‘\ IdT dt Then, d 0 2 dadt ( — g) (dT (dT (dT N + N gj da da dt ) — - da da’ 1 da — Vdt dt \. dt dt) Using the quotient rule for derivatives: dT da V di \. dadt) di dt \dt dadt di di / ,12 2 d = - da dt dadt dt dt dt dp dt) t’_(dadtdt).N (dt /4.4z The second term drops out because Therefore UI da \dadt N dl’ N•— dt — Vdt = di 1 d N..1 di Vdt dt but Using quotient rule N = 0 119 Appendix C. Derivation of Rate of Rotation of Generator Differential Equation — dTN dt — V dt IV di ) 2 dt di \. dt (dd \dt dt — (.& 2 \di di) 2 d dt 2 Vdt (E.& \di di) dt dt2 di) dt / dt Again, the second term drops out because dp dt — N=0 (ç.N) dTN dt Vdt di (. dadi g dT da (.& dadi dt) — /2.42 dt (42 Vdt ‘di / dp at di) (. \dadi dT a di Vdi Combining terms yields, • g) — -‘-.N i..dt2 ‘dadi J 1 d0 2 dadt ( \dadi g (42.2’ \di di) Using Identity: (A. C)(B. D) = (A x B) (C x D) yields, . ( d 0 2 dadt X • (N x g) (2 \dt di (x)gxN \cit — (dadt di ) 2 dt dt \dt But, A•(BxC)=(CxA).B — ) 2 ‘dt \di — — 2 \.di dadt di di 3 (& di di) d p 2 dadt 2 di (1 di di — (A. D)(B. C) 120 Appendix C. Derivation of Rate of Rotation of Generator Differential Equation 9 2 d — dt2 dt) . — dt dt) 121 (C.6) Appendix D Tension Catmull-Rom Spline The Catmull-Rom Spline matrix with a tension parameter, /3, is shown below. The phantom point end conditions are shown on the following pages. —2.0/3 4.0 1 P(t) = 3 u u 2 u 1 — 4.0/3 2.0/3 2.0/3 2.0/3 — 6.0 6.0 — — 4.0 2.0/3 4.0/3 —2.0/3 1 —2.03 0.0 2.0/3 0.0 0.0 2.0 0.0 0.0 Pi-’ Pi 1 Pi+ /3 = Tension parameter 122 Appendix D. Tension Catmull-Rom Spline D.1 123 Phantom point, P(O), at 0 —2.0/3 4.0 4.0/3 2.0/3 1 P(t) = — 2.0/3 2.0/3 — 6.0 6.0 — — 4.0 2.0/3 4.0/3 —2.0/3 0.0 0.0 0.0 1.0 —2.0/3 0.0 2.0/3 0.0 0.0 2.0 0.0 0.0 1 PiPi = Fl_i Fi+i 2 Fl+ 1 P 1 P 0.0 1.0 0.0 0.0 = Pi+i 2 Fl+ Fl__i Appendix D. Tension Catmull-Rom Spline —2.0/3 4.0 1 F(0) = 0.0 0.0 0.0 124 2.0/3 2.03 — 4.0/3 2.0/3 — 6.0 6.0 — — 4.0 2.0/3 4.0/3 —2.0/3 1.0 —2.0/3 0.0 2.0/3 0.0 0.0 2.0 0.0 0.0 Fi Pi 1 Fi+ i+2 Appendix D. Tension Catmull-Rom Spline D.2 125 Phantom point, P(1), at n —2.0/3 4.0 4.0/3 2.0/3 1 P(t) = — 2.0/3 2.0/3 — 6.0 6.0 — — 4.0 2.0/3 4.0/3 —2.0/3 1.0 1.0 1.0 1.0 —2.0/3 0.0 2.0/3 0.0 0.0 2.0 0.0 0.0 Pi-’ Pi 1 Pi+ 2 Pi+ Pi-’ Pi 0.0 0.0 1.0 0.0 = 1 Pi+ 2 Fi+ Appendix D. Tension Catmull-Rom Spline —2.0/3 4.0 1 P(1) = 126 — 4.0/3 2.0/3 2.0/9 2.03 — 6.0 6.0 — — 4.0 2.09 4.0/9 —2.0/3 1.0 1.0 1.0 1.0 —2.0/3 0.0 2.0/3 0.0 0.0 2.0 0.0 0.0 1 Pi1 P 1 Pi+ 1 Pi+ Appendix E Beta-Spline P(t) = [ 2 1 1 ] —2.0/3? 2.0(132 +/? +i3? + /9i) —2.0(132 +i3? +/3i + 1.0) 2.0 6.0/3? —3.0(132 + 2.0/3? + 2.0/3?) 3.0(132 + 2.0/3?) 0.0 6.0/3k 0.0 (5 —6.0/3? 2.0/3? 6.0(/3? (/32 — + 4.0/3? + 4.O/3) 1 P 1 P 1 Pi+ = 132 +2.0/3?+4.0/3?+4.0/3i+2.0 = /32 Bias Tension 127 2.0 0.0 Appendix E. Beta-Spline E.1 128 Phantom point, P(O), at 0 P(t) —2.O/3 = [0.0 0.0 0.0 2.0(/32 + /3 + i? + i3i) 1.0] —2.0(,62 + /? + /3 + 1.0) 2.0 6.0/9 —3.0(/92 + 2.09 + 2.0/3?) 3.0(132 + 2.0/3?) 0.0 /3) 6.03 0.0 + 4.0/3? + 4.0/3k) 2.0 0.0 S —6.0/3? 2.0,6? 6.0(13? (/32 — 1 PiPi 1 =Pi1 Pi+ 2 Fi+ S = 2 / ? 3 +4.0,B?+4.0i91+2.0 +2.0/ 40 Appendix E. Beta-Spline 129 1 PiPi 1.0 g 2.0/3 (/32 + 4.0/3 + 4.0/3k) 2.0 0.0 = i+1 2 Pi+ — (/32 )P + 2.0P÷ 1 1 + 4.O/9 + 4.O/3 (6—2.O,8) Appendix E. Beta-Spline 130 P(0) —2.O/3 1 S = [0.0 0.0 0.0 1.0] 2.0(132 + i? + I? + /3k) —2.0C82 + /? + th + 1.0) 2.0 6.0/3? —3.0(132 + 2.0/3? + 2.0/3?) 3.0(82 + 2.0/3?) 0.0 /3) 6.0,3 0.0 + 4.0/3? + 4.Ofli) 2.0 0.0 —6.0,3? 2.0/3? 6.0(13? (/32 — +2.oP )P i (I32-I-4.OI3?+4.o31 1 S—2.O3? Pi 1 Pi+ Appendix E. Beta-Spline E.2 131 Phantom point, P(1), at n P(t) —2.0,6 1 S = 2.0(/32 + /3 + [1.0 1.0 1.0 1.0] /? + 3) —2.0(132 + 6.0/3 —3.0(132 + 2.0/3 + 2.0/3?) —6.0/3 2.03 6.0(/3 (/32 /? + i3 + 1.0) 3.0(/32 + 2.0/3?) 0.0 /3i) 6.0/3k 0.0 + 4.0/3? + 4.0/3k) 2.0 0.0 — 1 PiPi = 2 Pi+ 1 Fi+ 2 Fi+ S 2.0 = /92+2.0/3+4.0/3?+4.0th+2.0 Appendix E. Beta-Spline 132 Pi-’ Pi 1.0 0.0 2.0/3 (4.0/3? + 4.0/3k + /32) 2.0 = 1 Fi+ i+2 — i+2 2.0/3F, + (4.0/3? + 4.06 + 11 )P 2 f3 8—2.0 2 Pi+ Appendix E. Beta-Spline 133 P(1) —2.0/3? 2.0(/32 = [1.0 1.0 1.0 1.0] + i? + i? + /3) —2.0(/32 + /? + /3 + 1.0) 2.0 6.0/3? —3.0(/32 + 2.0/3? + 2.0/3?) 3.0(/32 + 2.0/3?) 0.0 /3k) 6.0/3k 0.0 + 4.0/3? + 4.03) 2.0 0.0 45 —6.0/3? 2.0/3? 6.0(/3? (/32 — 1 Pi1 P 1 Fi+ 11 2.Of3?1’ + (4.O+4.O31+P2 )P 6—2.0 Appendix F Derivation of Normal Directrix Control Vertices It was found necessary to solve for a desired number of normal directrix control vertices from two space curves. The two space curves must be the same type of splines as the normal directrix in order to minimize the complexity of calculations. The number of control vertices of each space curve can be different from each other as well as both numbers can be different than the desired number of normal directrix control vertices. In order to solve for the normal directrix control vertices we solve a relation involving p(t), our known spline, r(t), the spline we wish to solve for, in the following error equation: 1 N -1 = error Ip(t) — 2 dt r(t)1 (F.1) where we wish to minimize the integral and solve for the Contro Vertices, C: JN1 öerror = 2 N—2 = = J’ dp(j+ s) - ld (p(t) — r(t)) dt (F.2) (+s) (p(j+s)—r(j+s))dS dC, 0 (F.3) (F.4) dp(j+s) pj + s)dS 2 j ’ = 134 . r(j + s)dS (F.5) Appendix F. Derivation of Normal Directrix Control Vertices 135 where the the above terms are defined below in vector form: Ci-’ Ci p(j + s) 3 = 2 [A] . (F.6) 1 Ci+ 2 Ci+ ij—1 6 [ ãp(j+s) = 2 1] [A] (F.7) Where A refers to the 4 x 4 spline basis matrix and 6 is the Dirac delta function. Using the identity (BA)T CT (AB)T = = ATBT expanding for triple products (ABC)T = CT BT AT we get the following relation: T 2 C ãerror [A]T [3 dS (F.8) 2 = s 1 2 C Appendix F. Derivation of Normal Directrix Control Vertices 136 1 N—2 = f r(j + s) [3 2 0 ] s 1 [A] dS (F.9) ij+1 5 ‘ ij+2 j Simplifying some of the terms yields the following form F j_’l 1 N2 ôerror ac T = I [A]f I ij+l 5 [3 2 S 1 I I ] dS[A] I (F.lO) [ [1 Sij+2j 1 Co 1 I 3=01 r I CN1] T ijl 5 ‘ — iii N_i I = 1 L 76541 I j=0 T [A] I . 1 43I [A] 1 1 1 1 1 5432 1j+2 6 1 1 432 r 0 Ic II I II I 1 I II (F.11) j N—2 (i,j)1’ = j=0 (F.12) Appendix F. Derivation of Normal Directrix Control Vertices 137 Co where 1’ (F.13) CN1 N—2 (i,j) [1’] = (F.14) j=o ij—1 5 N—2 1 r(j+s) dS s 1] [Al (F.15) ij+1 5 ij+2 8 (F.16) —1 3 i 5 ‘ N—2 [Fj = —1 [(ii)] N—2 1 j r(j+s) [3 2 dS (F.17) 1] [Al 1 sj+ ij+2 6 Once the control vertices have been solved for, the desired end condition constraints must be invoked. Namely: Co = Po+(Ci—Fo)g*g (F.18) CN = PN+(CN1—PN).g*g (F.19) Appendix F. Derivation of Normal Directrix Control Vertices 138 where P are space curve control vertices and C are solved for control vertices. Also, (P F 1 ) 0 1 0 IF’—P - = (FM gN = — PN) FM—FoI (F.20) (F.21) Appendix G Modified Conventional Approach Derivation A few trials of the normal directrix method yielded results which were very dependent upon the initial position of the normal directrix when trying to match the normal directrix method developable surface to two space curves. Further testing revealed that if a good initial guess of a normal directrix could be achieved to match to two space curves, then little optimization or correction is necessary and a closest developable surface could be found. Below is a figure relating the variables used in the derivation of solving for an initial normal directrix control vertices given the control vertices of two space curves. I.t)g(v) Figure G.1: Orientation of Space Curves and Directrix 139 Appendix G. Modified Conventional Approach Derivation 140 From Figure G. 1 f(u) and g(v) refers to the space curves and p(t) refers to the normal directrix. We relate these as follows: p(t) = (1—a)f(u)+ag(v) (G.1) Using partial fraction expansion from equation G. 1 we get: = (g—f) (G.2) = dfdu dgdv .da (l_a)T + a--- + gen (G.3) From Figure G. 1 we also calculate the normals at these points, namely: nl = ---Xgen (G.4) = dg j—Xgen (G.5) an Now, relating the out-of-plane curvature with the in-plane curvature of the two space curves we get the following: -; X 1 n ) . — du dt = I 211 (gen x 2 n ) . Idv (G.6) — dvi Also, relating the number of control vertices between the three curves and how and are related we have: 1.0 du 1.0 dv 2n dt +- 1.0 (G.7) Rearranging, gives: dv dt — — 2n 10 — nt du 2n dt (G.8) Appendix G. Modified Conventional Approach Derivation Rearranging, Equation G.6 to solve for and and f 2 d —1 II i• t gives: . ____ du 141 1 nt (genxni). II S g 2 d (G.9) II fl2 -II dg I (genxn ) 2 .— ) dv’ / dv — 2n,, ndu - da dt (G.1O) df dg d (1—a) . 1 gen+aj---.gen (G.11) Appendix H Relating the Space Curves to the Developable Surface Two forms were developed: 1. The first sets up the relations and solves for the equations 2. The second calculates the parameters directly The First Approach There are primarily three constraints which determine the distance from a point on the developable surface to the closest point on a space curve. They are as follows: dp(t) n(t) x g(t) q(t, s) = p(t) + sg(t) r(u) = q(t, s) + ln(t) where, n(t) is the normal of the directrix, q(t,s) is a determinable point on the surface, resulting in the determination of the corresponding location along the space curve. There are three variables in which we need to solve for for a corresponding value of t, which is the key parameterized variable. In order to solve for these we try to solve the problem in terms of a differential equation relating the variables u,s,l, to t. This is done 42 Appendix H. Relating the Space Curves to the Developable Surface 143 by partial differentiation of the function r(u). The differential equation is with respect to t is as follows: drdu dp dg ds dn dl = drdu Rearranging into a useful form, dp dg dn ds dl -g--n = There are only three unknowns to solve for given the initial conditions relating each space curve to the developable surface: 1) At t = 0.0, 11 2) At t = 0.0, ul 3) At t = 0.0, q(0.0,sl) 4) si = s2 (q(0.0, s2) = = = (q(0.0, si) — 0.0 and 12 — = 0.0 and u2 0.0 0.0 = rl(0.0) and q(0.O,s2) r2(0.0) = p1(0.0)) g(t) p2(0.0)) g(t) and, and can be evaluated, dg dt dn - - — — (Y)dp \dt p 2 d 1 dtJ dp 1 dg Xgj+X Three unkowns are then solved for each space curve and the developable surface. Appendix H. Relating the Space Curves to the Developable Surface F dux —g X du dt L duy —gy fly ds dt duz —g nz dl dt I I R dty 144 +1411 +s dty dty dL x The Second Approach This approach uses the same definition but only differs after the follow differential equation has been derived. drdu ds dl - g- n — (11.1) Taking the dot product with vector g: ( du — ds g g dl — n g (11.2) dt g +s- dn g +l- This simplifies to the following: ds -ï (H.3) (11.4) ( ...I’du (11.5) Similarly, dot product with n, (H.6) and simplifying yields: (H.7) ( dl j’\du (11.8) = Similarly, (11.9) dot product withy (11.10) ( 1\ dp’du dt) dt ;; +ldt _.; — dt dt+Sdt dt dt (H.11) Appendix H. Relating the Space Curves to the Developable Surface 145 Note perpendicular (11.12) vector relationships (H.13) Forum = ;_ dt — dt (11.14) g (11.15) g 2 dt (11.16) — / —.. .12 (;.; 1 dt but, identity:(A x dt = xg (CxA)B B) C —.-. (dp / therefore, x dt therefore, — (c x = x _; Y) ( ;dt du ) yielding: (H.20) _.g+l x g (H.21) IT i•_ (11.22) 2 dt — - — ( 2 d1 (11.23) — dp dp dt dt / du (11.19) I x g . (H.17) (H.18) (BXC)A B) C dt —s. dp\dg dt = ( & but Identity:(A = _—.—.\ iç Simphfyrng: dt +x — .—.—. dp dp ..12 2 dt P — I • g —‘-• fl (11.24) (11.25) p 2 d ._.(sg = +ln)) (11.26) (dr dp’\ Appendix I Intersection of Developable Surfaces Once the developable surfaces have been created there is then the requirement that they must intersect without gaps or spaces between them. Shown below is the vector representation of the initial conditions and their orientation relative to each other. In Figure 1.1 we see the various vectors in which their dependence to each other will be shown as follows: 4U t dr 1 dt V Figure 1.1: Intersection of Three Developable Surfaces 146 Appendix I. Intersection of Developable Surfaces 147 As shown in Figure 1.1 the independent parametric variable is t. The two dependent variables are u and v. There are three developable surface control vertices with all func tions of the independent variable t. The independent space curve is p(t); one dependent space curve which is a function of u is f(u); the other dependent space curve is g(v), which is a function of v. The intersection curves joining two developable surfaces are r 1 which is relating f(u) and p (t) and the other intersection curve r 2 which relates g (v) and p (t). Derivation of the equations used to calculate the intersection of the developable sur faces are as follows: = 1 dr 1 dr (df ---_xg = = p+si: dp - + si (1.1) dg -- + dsit -- g (1.2) 0 (1.3) (1.4) = dp (df L\ 1 ds dt — - — — du 1 r = gj / dg +s1-- (df X (15) d p 2 . dt d -— dt 16 (.) !. (17 di f 2 d uju w dt du ugu 2 f+s du (1.8) Appendix I. Intersection of Developable Surfaces 148 dridridu dr 1 du dtdt dr df — du (1.9) + dg dsu) + g— du du du (dhdfu + u — du Since (Ho) S2u— dg S2u du — + ds2U du is perpendicular to du the last term drops out Substituting gives, _- ( — 1 dr du df df • — df — du (1.11) (1.12) (1.13) — ‘‘ df S2u (1.14) = But, (1.15) 1 r = 2ugu = and, r 1 = 4 + s2ugu, and f, fu (1.17) + s1tg (1.18) — p (1.16) Substituting gives, (1.19) f 2 d - u du dffudu(dhdh du f 2 d _.(ri_fu)) = (1.20) (1.21) du(dh4fu = d2f — _.T_T.(P+s1_fU)) (1.22) du(dh du (; dgdgdu ----.i = P + (1.23) (1.24) (1.25) Appendix I. Intersection of Developable Surfaces dr 2 dt 1 dr (dhv — — ; 149 2 ds (1.26) =0 I\d — — (1.27) dp dt dg (dh idv +S2t dp Idh (dh dg — -;-. 2 +S •. ds i 2 dt Idh I\d — x + ( x (1.28) - (1.29) — X g, p 2 d •iti•Yt — dg dt — dg dt — — dp dt dp dp dt dt h 2 d —-•g (1.30) — 2 r (1.31) - = dr 2 dt dr 2 dt dh dv dv h+si (1.32) dr dv 2 dv dt dr 2 dt = dh dv = dv (1.33) I dh + dv (dh Since dg = Substituting for 2 r = S1,j9.j = r2 = h+si + g, dsi dg dh . dv s dh (1.34) 2 ds 0 the last term drops out gives, dh (1.35) (1.36) (1.37) (1.38) (1.39) p+S2tTh (1.40) Appendix I. Intersection of Developable Surfaces dr2dh = . -— . 150 (1.41) - dv(dhvdI d2h = (r2 — — )) 0 h (1.42) dv(dhvd1 d2h = ( (1.43) dv(dh0d’ dg + -- dv + — dgdgdv - . (p + (1.44) 2 S (1.45) Appendix J Derivation of Flat Plates Once a developable surface has been created, the shape and position of the consecutive flat plates that make up the developable surface need to be known. Below is a simple derivation of the algorithms used to provide the information of the shape of the plates in the x-y plane and the angle and rate of change of bending of the plates relative to each previous one. In order to show the relations we need to define the various vectors in the two different spaces. In the first space we are in a three dimensional system with the following termns: The vector tangent at a point along the directrix is defined as at a point along the directrix is defined as 3 dimensional directrix is defined as . . The vector curvature The vector generator at a point along the 93d. In the second space we are in a two dimensional system with the following terms: The vector tangent at a point along the directrix is defined as ture, to be integrated to find the next tangent, is defined as a point along the 2 dimensional directrix is defined as . . The vector curva The vector generator at g2d. Since we are resolving the desired information in the 2 dimensional x-y plane, the normal, n, is obviously n = (0,0, 1). Setting up the differential equation we intend to solve, we resolve the vector relation ships into either in-plane or out-of-plane components. The following differential equation 151 Appendix J. Derivation of Flat Plates 152 uses the existing definition of the differential equation defining the calculation of a gen erator and relates that as the in-plane component which defines the first term of the differential equation. The out-of-plane component is just resolving the out of plane cur vature of the directrix and the three dimensional generator from the three dimensional system. The resulting differential equation is used to find the new tangent and corresponding point on the x-y plane where the next point of the plate is located: 2 dq = 2 ‘d± dt)dt () /2 dp +-g3d g2d Appendix K O.D.E. Class, Adaptive Step Size and Runge-Kutta Method Previous attempts at using different numerical integration methods proved that a fair amount of processing overhead takes place when making a function call to a numerical integration subroutine. Arguments must be pushed onto the stack, a call is then executed, a return is then implemented concluded with the parameters of the stack being removed. Instead of this traditional approach, the numerical integration function was placed “inline”. This technique is used in such languages as C++ and declared as an “inline” function. One may argue that this is sacrificing “modularity”. It is, but when more than one function call is being made, like in a loop, then using an “inline” function shows considerable gains in speed, smaller executable files, and less functional overhead. Evaluating the Initial Value Problem (I.V.P.): y’=f(x,y) where y(x )=y 0 Runge-Kutta formula involves a weighted average of values of f(x,y) taken at different points in the interval: Xfl X XnI1 The fourth-order Runge-Kutta is written in one of the classical forms as follows: Yn+i = y, + h )+m 3 ) 4 + 2(m2 + m 153 Appendix K. O.D.E. Class, Adaptive Step Size and Runge-Kutta Method where m 1 The sum = f(x,y) 2 m = f(x + 3 m = f(x + i-h, y, + hm ) 2 4 m = f(x + h,y + hm ) 3 +2m3+m 2 (ml+2m ) 4 154 + hmi) can be interpreted as an average slope. m 1 is the slope at 2 is the slope at the midpoint using the Euler formula to the left-hand of the interval, m go from x to x + , m is a second approximation to the slope at the midpoint. Finally, 4 is the slope at m Xr, 3 to go from x to + h using the Euler formula and the slope m x + h. This is a very accurate formula ( halving the step size reduces the local formula error by the factor it is not necessary to compute any partial derivatives of Another note should be mentioned, that if f does pj depend on y, then: 1 m = f(x) 2 m = 3 m 4 m = f(x + h) = f. 1 f(x + h) and the equation reduces to that of Simpson’s rule when evaluating the integral of y’ =f(x): Yn+1 — Yn = [f(n) + 4f(x + h) + f(x + h)] 5 which is in agreement with error in Simpson’s rule has an error proportional to h Runge-Kutta formula [6]. Appendix L Non-Linear Optimization Downhill Simplex Method The Downhill Simplex Method is a non-linear optimization tool for multidimensional min imization problems, ie. finding the minimum of a function of more than one independent variable. The original paper citing this method can be found by the authors Nelder and Mead [16]. A more recent overview with code samples can be found in Press [20]. A quick overview of the method and code implemented is presented here and mentioned briefly in Section 4.3.1 in this thesis. A simplex is the geometrical figure consisting, in N dimensions, of N + 1 points (or vertices) and all their interconnecting line segments, polygonal faces, etc. In three dimensions it is a testrahedron, not necessarily the regular tetrahedron. The algorithm is supposed to then make its own way downhill through the geometrical configuration of an N - dimensional topography, until it encounters an (at least local) minimum [20]. The downhill simplex method must be started with N + 1 points defining an initial simplex. It then starts and takes a series of reflections, contractions and expansions. Each one of these is assigned a variable: cv, reflection; /3, contraction; y, expansion. The first variable, a, is a positive constant, which is a scalar multiplication constant which mirrors the point through the simplex. The second variable, /3, is a constant which values lie between 0 and 1. It is a ratio of the distance of the point relative to the simplex centroid. The third variable, y, is greater than unity and it is a ratio of the current point to the centroid with a point along the line joining the point to the centroid. In figure L.1 illustration a) shows the variable a being implemented in the algorithm. 155 Appendix L. Non-Linear Optimization Downhill Simplex Method 156 Illustration b) shows both reflection and expansion. Illustration c) shows a contraction. Downhill Simplex Method (a) (b) (c) Figure L.1: Analogy of a Simplex for the Downhill Simplex Method A few trial runs were done varying the three variables, c, table L.1 below. *define ALPHA 0.9 *define BETA 0.4 *define GA11IA 1.9 *define GET_PSUW for (j0;j < ndi.;j++) su + pLi] [j]; { for (i0,sii.0.0;i<mpts;i++)\ psu[j]su;} , -y and are listed in the Appendix L. Non-Linear Optimization Downhill Simplex Method a/37 2 4 2 3 1 3 20 1.0 4.0 1.0 2.0 0.4 0.9 1.9 Table L.1: Values of variables used for Simplex Method void si.plex(double **p,double *y,iat ndia,double ftol,double (efunk) (double *), hit *nfunk, hit nan) { mt mt i,j,k,ilo,ihi,hthi,apts—ndiatl,rnnout; coiuttO; double ytry,ysave ,sua,rtol,*psua; pent = new double[ndia]; for(i0; i < ndia; i++) psna[EfrO.O; double aaotry(double ** ,double * ,double * ,int,double (*)(double *), mnt,int *,double); *nfunkO; GET_PSUM { for (;;) iloO; ihi = y[O]>y[1J ? (inhil,O) for (i0;i<apts;i++) if (y[i] < { ylio]) ioi; if (yti) > ytihi)) { inbiihi; ihii; } else if (y[i] > y[hthi)) (inhiO,1); 157 Appendix L. Non-Linear Optimization Downhill Simplex Method if (1 ihi) inhii; } rtol2.0*fabs(y[ilui)—y[ilo])/(fabs(y[ihi])+fabs(y[ilo))); if (rtol < ftol){ break; } if (enfunk > unx){ coutW’\nToo zany iterations in SINPLEI\n”; goto runout; } ytrya.otry(p ,y ,psu.,ndia,funk,ihi,nfunk,-ALPRA); if(ytry C 0.001) tol = 0.0000000000001; cout<<”\nl if( kbhitO) goto runout; if (ytry <S yEilo]){ ytrya.otry (p ,y ,psu.,ndia,funk, liii ,nlunk ,GMIKO; if(ytry C — .001) tol = 0.0000000000001; cout<<”\n2: “<<ihi<<’\n”; if( kbhit() ) goto runout; } else if (ytry > y[inlti]) { ysavey[ihi]; ytryamotry(p,y,psia,ndi.,funk,thi,nfunk,BETA); if(ytry < 0.001) tol = 0.0000000000001; cout<<”\n3: “<<ihi<<”\n”; if( kbhit() ) goto runout; if (ytry > ysave) { for (i0;iCapts;i++) if (i ilo) { { for (j0;j<ndi.;j++){ psu.Ej)0.S*(p[i] fj]+p[ilo] ph) [j]psla[j); } y[iJ(*funk) (psusO; if(kbhitO) goto runout; } } *nfunk + GET_PSUM ndia; [JU; 158 Appendix L. Non-Linear Optimization Downhill Simplex Method } } :i runout: rmtouto; delete pens; } double asotry(double **p,double *y,double tpsns,int ndia, double (*funk)(double *),jnt thi,int *nfuDk,double fac) { Ant j; double fad ,fac2,ytry,*ptry; ptry = new double[ndiml; facl(1 .O-fac)/ndi.; fac2facl-fac; for (j0; j<ndia; j++) ptry[j]psus[j) *facl—p[ihi] fj)*fac2; ytry(sfinik) (ptry); if (ytry < yfihi]) { yEthi]ytry; for (j0;j<ndia;j++) psua[j1 += { ptry[j]—p[ihi] [j); p[ihil [jfrptry[jl; } } delete ptry; return ytry; } tundef ALPHA *undef BETA tundef GAMMA 159 Appendix M Codelisting M.1 Class Tools Used M.1.1 vector.h *tnclude <iostreea.h> *ifndef _VECTOR_ *define _VECTOR_ enun direction { I, Y, Z }; class zatrix; class vector { lat n; static mt default.lengtb; float selenent; public: vector() {ndefault_length; elementnew float [n) ;}; vector(int din) {ndia;elenentnew float[n];}; vector(int din, float x); vector( coast vectort v ); vectorQ; static void set_default(int n){default_lengthn;}; static void set..AefaultO{default_length3; }; float& operatorEl( tnt i) {return elesent[i];}; float operator[](mnt i) coast {retnrn eleaent[i] ;}; friend tnt size(const vectork v) { return v.n;}; vectort operator(float x); vectort operator(const vectork v); 160 Appendix M. Godelisting 161 vectork operator(const matrixk a); friend vector operator+(const vectort a, const vectort b); vectort operator+( const vectort a); friend vector operator—(conat vectort a, coast vectort b); vectort operator—(coust vectort a); friend vector operators(const vectort v, float a); friend vector operators(float a, const vectort v); friend float operators(const vectork a, const vectort b); II dot product operator friend vector operators(const matrixi a, coast vectort b); vectort operator*(float a); vectort operatorn(const vectort a); II elaent—by—eleaent multiplication friend vector operator/(const vectort v, float a); friend istreant operator>>(istreaak, vectort); friend ostresat operator<<(ostresat, coast vector&); tendif M.1.2 matrix.h tinclude <vector .h> #include <iostream.h> tifndef _NATRIX_ *define _MATRII_ class matrix mt { nr, nc; float **element; public: matrix(int rows, mt columns ); matrix(const matrixi a); matrix(const vectort a); ‘matrixO; matrixt operator(float x); matrixi operator(const matrixi); float* operatorO (hit i) {return element[i] ;}; coast floats operator [1 (tnt i) coast {retuxa element Li);); vector row(int i); friend hit n_rows(const matrixi a){return a.nr;}; Appendix M. Codelisting friend mt 162 n_coluns(const matrix* a){return a.nc;}; vector column(int i); friend matrix exp(const aatrixk a); friend matrix elaent_ault(conat matrixi a, const .atrixt b); friend vector diagonal(const matrixt a); friend matrix operator+(conat matrix& a, const matrix a); friend matrix operator—(const matrixk a, cost matrixt b); friend matrix operators(const matrixk a, cost matrixi b); friend vector operator*(const matrixk a, cost vectort b); friend matrix operator*(const matrix& a, float b); friend matrix operator*(float a, cost matrixi b) {return b*a;}; friend matrix transpoae(conat matrixt a); friend matrix solve(const matrix& a, cost matrix& b, float edet friend matrix identity(int); friend matrix inverse(const matrixk a, floats detJULL); friend istream& operator>>(istream&, matrix&); friend ostreaak operator<<(ostream&, cost matrixi); friend mt operator<(cost matrixi, cost matrixi); friend matrix abs(const matrix&); #endif M.1,3 develop.h tinclnde ‘matrix .h” tinclude <fstream hpp> *ifndef _CURVE. #def Inc _CURVE_ class Curve{ hit n; char splinetype; vector *point; static matrix BB; public: void aptype (char what) {splinetypenhat ; }; char typeofsplineO{return splinetype ;}; NULL); Appendix M. Codelisting 163 tnt sizeO{return n;}; CurveO{n5; splinetype ‘1’; point = new vector[5];}; Curve(int nua){n=nus; splinetype’l’; point = new vector[n];}; Curve(Curve& a); CurveO{delete[n] point;}; void resloc(int nun); vectork operator[] (mt i){return point[iJ ;}; Curvet operator=(Curve& c); vector operator() (double t); vector tsngent(double t); vector curvature(double t); double deriv_2c(double t, double deriv_c(double t, mt i); mt i); static void initbasis(char sptype’2’, double b11.0, double b20.0); static aatrir EBB; tendif *ifndef _SURF_ tdefine _SURF_ class Surface{ private: Curve eleftedge, erightedge; Curve *directrix; Surface *leftsurf, *rightsurf; public: Surface( Curve *leftedge, Curve srightedge, Surface eleftsurf, Surface *rightsurf ,int ncontrolpnts5 ,double step&T .0, double outpo.05,int surfnot); void Rightsiuit(Surface trights); void Leftsinit(Surface Clefts); double Opttnize(double toleranceo.2, tnt iter.ar = 250); void Develop(void); void Developt (void); void Developtl(void); void DraIL3d(int genneso,double stepl.0,mnt surfnol); void Draw_33d(int gennoso,double stepl .0, tnt surfnol); void Draw_3id(int genno5o,double step7.0,tnt surfnei); Appendix M. Codelisting void Draw_2d(iat genno=5O,double steop7.O,int surfao=1); surfaceO; Surface(Surface& s); Surfacet operator(Surface& s); Sendif tifudef _FILELISTS_ Idefhte _FILELISTt tinclude <fstreaa.hpp> class Filelisto{ public: ofstreaa file; Filelisto *nezt; ofstreaat operatorO(int tO; class Filelisti{ public: if stream file; Filelisti tusit; ifstreamt operatorlj(int n); tendif M.1.4 ode.h *ifndef _ODE_ *define _ODE_ tinclude “vector .h” class dynamic_systa { private: double time, step_size; vector state; vector error_scale; vector (ederivative) (double, vectort); 164 Appendix M. Codelisting 165 public: dynaic_syste.(vector (*)(double, vectort), double start_time, vectort initial_state, vector *error_scaleo); doublet whenfl; void resetO; doublet operatorD(int i); vector operator() (double t); vector step(double delta); to, void rk(double doublet del_t, double ti, vectort z, vector (*f)(double t, vectort x), vectort err_scale); tendif M.1.5 vector.cpp tinclude <ioatreaa.h> tinclude <process .h> S include <atrix h> const char SP = static inline mt ‘ sin(int a,iut b) { return a>b?a:b; } mt vector: :default_length = 3; vector::vector(iut dim, float initial_value) { ndim; eleaentnew float [it]; for(int i0; i<n; i++) elaent[i] = initial_value; } vector::vector(const vectort a) Appendix M. Codelisting { n = a.n; element = new float[n); for(int 1=0; i<n; 14+) element [1) a element [ii; } vector : :vector() { delete element; } vectort vector: : operator(float r) { for(float* te.pele.ent+n-t ;temp>element ; teep——) *tempx; return ethis; vectort vector::operator(const vectort v) { 1f(nn.n){ delete element; nn n; elementnew float [n]; } for(Int 1=0; i<n; j++) element [1] n element [1]; return ethjs; } vectort vector::operator(const aatrix& a) { if(n!n.rows(a)){ delete element; nn_rows (a); element = new float [n); I for(int 1=0; i<n; j++) 166 Appendix M. Codelisting elaent [i] = &• ele.ent fi] [0]; return *this; } vector operator+(const vectort a, const vectort b) { vector teap(a); for(int i0;i<.in(a.n,b.n) ;i++) tep.elesent[i] + b.eleaent[i]; return teap; } vector operator—(const vectort a, const vectort b) { vector temp(a); for(int i0;i<ain(a.n,b.n) ;i++) te.p.element[i] —= b.eleaent[i]; return teap; } vector operator*(const vector iv, float a) { vector b(v); for(int iv.n—1;i>0;i——) b.eleaent[i] *= a; return b; } vector operators(float a, const vectori v) { vector b(v); for(int in.n—1;i>0;i——) b.eleaent[i] s a; return b; } float operator*(const vectori a, const vectori b) { float tempo.O; 167 Appendix M. Codelisting for(int i’ain(a.n,b.n)—l;i>O;i——) temp + a.eleaent[i]sb.ele.ent[i]; return teap; } vector operator/(const vectort v, float a) { vector b(v); for(int in.n—1;i>0;i——) b[i] / a; return b; } vectort vector::operator+(const vectort a) { for(int 1min(n,a.n)—1;i>0;i——) element[i] + a.element[i]; return *thjs; } vectort vector::operator—(const vectort a) { for(int izmin(n,a.n)—1;i>O;i——) element [1] —= a. element Li.]; return *thjs; } vectort vector: :operatorn(float a) { for(int in—1;i>0;i——) ele.ent[i] *= a; return ethis; } vectork vector::operatorn(conzt vectort a) { for(int imin(n,a.n)—1;i>=O;i——) ele.ent[i] fl a.ele.ent[i]; 168 Appendix M. Codelisting return ethis; } istreamk operator>>(istrea.& input, vectort a) { for(int j0;i<a.n;i++) input>>a clement Li]; return input; } ostream& operator<<(ostresmk output, conat vectort a) { for(int i0;i<a.n—1;i++) output << a.element[i] << SP; output << a.element[a.n—i1; return output; } M.1.6 matrix.cpp #include <matrix .h> Siuclude <iostream.h> tinclude <process .h) matrix::matrix(int rows,int columns) { nccolumns; nrrows; elementnew float * Ens]; for(int 1=0; i<nr;i++) elesentli] = new floatlnc]; } matrix::matrix(const matrix& a) { mt i,j; nc = a.nc; nr = a.nr; 169 Appendix M. Godelisting element = new float * [nr]; for(i0; i<nr; i++){ element[i] = new float[nc]; j++) for(j0; j<nc; element [ii fjfra[il U]; } } matrix::matrix(const vectort a) { mt i; nc = 1; nr = size(a); element = new float *[nr]; for(i0; i<nr; j++){ element[i] = element[i][O] new floattnc]; = a[i]; } } matrix:: matrix() { for(int i0;i<nr;i++) delete element[i]; delete element; } matrixt matrix::operator(float x) { for(float** te.ptelement+nr—1 ; templ>element ; tempt——) for(float* temp2*tempt+nc—1 ;temp2>fltempl;teap2”) *temp2x; return *this; } matrixt matrix::operator(const matrixk a) { for(int i0;i<nr;i++){ for(int j0;j<nc;j++) 170 Appendix M. Codelisting clement [ii [ii = a. element [i] 171 [ii; } return *this; } vector matrix: :row(int i) { vector out(nc); for(int j0; j<nc; out[j] = j+) element[i][j1; return out; vector matrix: :column(int 1) { vector out(nr); for(int j0;j<nr;j++) outiji = elementEjiti); return out; } vector diagonal(const matrix La) { vector out(a.nr<a.ncra.nr:a.nc); for(int 1=0; i<size(out); i++) out[i) = a[i] Li); return out; } matrix element_mult(const matrix La, const matrixt b) { matrix result (a.nr,a.nc); for(tht i0; i<a.nr; i++) for(int j0;j<a.nc;j++) result Ci] Cj] = a[i] [j]*b[i] Li]; return result; matrix operator+(const matrixt a, const matrixk b) { Appendix M. Codelisting 172 aatriz te.p(a.nr,a.nc); for(int 10;i<a.nr;i++){ for(int j0;j<a.nc;j++) te.p ele.ent [1] . [j) = a. eleaent [1] [j]+b .elaent [1) [jI; } return te.p; } aatrix operator—(const •atrix& a, const .atrixk b) { .atrix teap(a.nr,a.nc); for(int 1=0; i<a.nr;i++){ for(int j0;j<a.nc;j++) te.p.element[i] [ii = a.element[iJ [jl—b.elesent[i1 [ii; } return teap; } aatrix operators(const aatrix& a, conat ntrixk b) { aatrix teap(a.nr,b.nc); for(int 10;i<a.nr;i++){ for(int j0;j<b.nc;j++){ float total = 0.0; for(int k0;k(a.nc;k++) total + a.eleaent[i) [k]*b.eleaent[kJ [ii; te.p.eleaent[i1 [jitotal; } } return temp; } vector operators(const matrix& a, const vectork b) { vector result (n_rows (a)); resulto.0; for(int i0;i<size(result) ;i++){ for(const float *tap_a& (a[i] [01), *teap_bb eleaent ; teap_b<b . ele.ent+b n ; te.p_a++,te.p_b++) . result [1] +fltap_a*tteap_b; . Appendix M. Codelisting 173 } return residt; } istreaa& operator>>(istreaak input, •atrix& a) { for(int i0;i<a.nr;i++){ for(int j0;j<a.nc;j++) input >> a[i][jl; } return input; } ostrea& operator<<(ostream& output, coast aatrix& a) { for(int i0;i<a.nr;i++){ for(int j0;j<a.nc—i ;j++) output << a[i][jJ << ‘ output << a[iJ[a.nc—i] << ‘\n’; } return output; } .atrix transpose(const ntrix& a) { satrix te.p(a.nc,a.nr); for(int i0; i<a.nc;i++){ for(int j0;j<a.nr;j++) temp.element[il fjfra.elementlj] } return te.p; } aatriz iclentity(int a) { satrix te.p(a, a); for(int i0;i<a; i++){ for(int j0;j<a;j++){ te.pEi] [j](ij)?i.O:O.O; Li]; Appendix M. Godelisting I return teap; I satrix operator*(const aatrixk a, float b) { tnt i,j; satrix te.p(a); for(i0; i<a.nr; j++){ for(j0; j<a.nc;j++) temp.element[il[j] *= b; I return teap; } tnt operatorfl(const aatrix* a, const ratrix* b) { tnt i,j; for(i0; i<a.nr;i++){ for(j0; j<a.nc;j++){ if(a[i] [j)>bfiJ [ii) retuni(0); I I return(i); } satrix abs(const aatrix& a) { .atrix teap(a); tnt i,j; for(i0; i<te.p.nr; i++){ for(j0; j<te.p.nc ; if(teap[i] teap[iJ[j) [ii <0.0) = I I return te.p; —teap[i][jJ; 174 Appendix M. Codelisting 175 } solve.cpp M.i.T S include <fstream . hpp> Sinclude <stdlib.h> tinclude “matrix .h” coast double TINY = i.Oe-35; inline void error(char * message) {cerr < inline double abs(double & a) message;exit(1);}; { return a<0.0?—a:a; matrix solve(matrix&a, matrix& b, double *det) { double *ptemp; if(a.nr!=a.nclla.ncgb.nr) error(”matrix::solve Attempted operation on incompatible matricies\n”); matrix c(a); double *scale; if(det !=IULL)*det=1 .0; scale = new double[a.nr); for(int i=0;i<a.nr;i++){ 1/loop over all rows, finding the largest double largestO.0; I/element in each row for implicit scaling. for(int j0;j(a.nr;j++){ double tempabs (c element [i] . [j]); if(temp>largest) largestte.p; } if(largesto .0) error(”matrix: :solve singular matrix\n”); scale[i]1 .0/largest; } for(int j0;j<a.nr;j++){ I/loop over all columns for(i0;i<j;i++){ I/solve for the elements of the upper triangular Appendix M. Codelisting double sum c.elementLi][j]; = 176 If matrix. for(int k0;k<i;k++) sun —= c . element Li] [k] *c . element I:k] c . element Li] [j] = U]; sum; } double largest mt 0.0; = imaxj; for(ij;i(a.nr;i++)( I/solve for the elements of the lover triangular double sum c.elementfi]Lj]; I/matrix. = for(int k0;k<j ;k++) sum —= c . element Li] [k]*c .element Uk] c element Li] Li] . = Li]; sum; double temp=scaleLi]*abs(sum); I/keep track of the largest element. if(temp>largest){ largest = temp; imax=i; } } if(imax!=j){ ptemp = f/if necessary, interchange rows. c.elementLj]; //IOTE: this row interchange method depends on c.elementLj] c.element[imsx]; f/the method of storing the matrix. = c.elementLimax] ptemp = = ptemp; b.elementLj]; b.elementLi] f/Interchanging rows of the RES matrix now b.elementLimax]; I/relieves the necessity of tracking the = b.elementLimax] = if(det!=NULL)*det ptemp; f/row interchanges for later use. = (*detl.0)?—l.0:l.0; scale Limax] ncale Li]; } if(c .elementLi] Li]0.0) c .element Li] Li]TIIY; double temp = 1.0/c.elementfi]Lj]; for(ij+i ; i<a .nr; i++) c.elementLi]Li] n temp; } f/The LU decomposition in now complete. if(det !IULL){ for(i0;i<a.nr;i++) f/Calculate the determinant of the matrix. *det } * c.elementLi]Li]; Appendix M. Codelisting 177 for(j0;j<b.nc;j++){ I/Solve for each colr of the b matrix. mt ii —1; = for(i0; i<a.ur;i++){ //Forwsrd substitution. double sum mt = b[i] [j); k; if(ii>—1) for (kii ;k<i ;k++) sum —= c.ele.ent[i] [k)eb[k) U]; else if(su.O.O) iii; b[i] [j]sum; } for(ia.nr—i;i>0;i——){ f/Back substitution. double sum = for(mnt ki+i;k<a.nr;k++) c.element[i][k]tb[k][j]; sum —= b[i] Li] sum/c.elaent[i] [i]; = } } delete scale; return b; } matrix inverse(matrixk a, double* det) { matrix b = solve(a, identity(a.nr), det); return b; } double det (matrixka) { double *ptemp; double deter = 1.0; matrix c(a); double Sscale; scale = new double[a.nr]; Appendix M. Codelisting 178 for(int i0;i<a.nr;i++){ i/loop over all rows, finding the largest double largesto.0; i/element in each row for implicit scaling. for(int j=O;j<a.nr;j++){ double te.pabs (c element Li) . U)); if (temp>largest) largesttemp; } if (large st 0 .0) return 0.0; scale Li) =1.0/largest; } for(int j=0;j<a.nr;j++){ //loop over all columns for(i0;i<j;i++){ i/solve for the elements of the upper trianguJLar double sum = c.elementLi)Lj); // matrix. for(int k0;k<i;k++) .elementLi) Lk]ec.elementLk) Li]; c.elementLi)[j) = sum; sum —= c } double largest mt 0.0; = i.marj; for(ij;i<a.nr;i++){ i/solve for the elements of the lower triangular double sum = c . element Li) U); //matrix. for(int k0;k<j ;k++) sum c . element Li) —= c.elementUi)Uj) = Lk) Cc. element Uk) U); sum; double tempscaleLi)eabs(sum); //keep track of the largest element. if(temp>largest){ largest = temp; imaxi; } } //if necessary, interchange rows. if(imaxj){ ptemp = c.elementLi); /iIOTE: this row interchange method depends on c.elementUj) = c.elementLimax) c.elementLimax); //the method of storing the matrix. = ptp; deter(deterl .0)?—1 .0:1.0; scale Lims-x)scale Li); } Appendix M. Codelisting 179 if(c .ele.ent[j] [j]0.O) c eleent [j] [j]TIIY; double temp 1.O/c.ele.ent[j][j]; = for(i=j+1 ;i<a.nr; j++) c.ele.ent[i][j] } *= temp; f/The LU decomposition in now complete. for(i=O;i<a.nr;i++) f/Calculate the deterninant of the deter natrix. c.element[i][i]; * return deter; } M.1.8 develop.cpp ///////////////,/////////////////////////////////////////////////////////////// // Class Curve C++ // // Modified Aug 8, 1992 by Brian Konesky, error noted deriv_2c tinclude <nath.h> tifdef __ZTC__ #include <fstrean.hpp> telse *include <fstrean.h> Sendif Sinclude “vector .h” finclude “natrix .h Sinclude “develop. satrix Curve: :BB(4,4); natrix Curve: :BBB(4,4); void Curve::realoc(int nu) { delete En] point; point new vector[nu); } Curve::Curve(Curve& a) { If Appendix M. Codelisting na splinetypea. splinetype; point= new vectortn]; for(int 1=0; i<n; i++) point[i]a.point[i]; } Curvet Curve::operator(Curvek c) { i:f(n ! c.n){ delete[n] point; nc.n; splinetypec splinetype; point = new vector[n]; } splinetypec. splinetype; for(int 1=0; i<n; i++) point[i]c[i); return *this; } vector Curve::operatorfl(double t) { mt t2; double t3 ,bbo ,bbl ,bb2 ,bb3; vector p; natrix pi(4,3),tt(1,4),c2(1,31,cc(1,4); bbO = BB [0] [0] +BB [1] [0] +BB[2] [0] +BB[3] [0]; bbl = BB[0] [1]+BB[1] [1]+BB[2] [1]+BB[3] [1]; bb2 = BB[0] [2]+BB[1] [2]+BB[2] [2]+BB[3] [2]; bb3 = BB[0] [3] +BB [1] [3]+BB[2] [3]+BB[3] [3]; t2(int)floor(t); t3t—(double)t2; if(t == n i){ — t3 = 1.0; t2 = n — 2; } tt [0] [0] t3*t3*t3; tt [0] [1]t3*t3; 180 Appendix M. Codelisting 181 tt[0] [2]t3; tt[0] [3]1.0; if(splinetype tf(t2 == == 0){ pi[O] [0)point[t2] [0]; pi[0] [1]point[t2] [ii; pi[0] [2]point[t2] [2); } { else pi[0] [O]point [t2—1] [01; pi[0] [1]point[t2—1] [1); pi[0] [2]point[t2—1] [2]; } if(t2 == a — 2){ pi[3] [0]izpoiat[t2+1] [0]; pi[3] [1]point[t2+1] [1]; pi[3] [2]point [t2+1] [2]; } { else pi[3] [0]’oint[t2+2] [0]; pi[3] [1]point[t2+2] [1]; pi[3] [2]point [t2+2] [2]; } } else if (splinetype if(t2 == == ‘2’) { 0){ pi[O][O] = ((poiat[t2][0])e(1.0 — BB[3][1]) —(point[t2+1] [0])*BB[3] [2])/BB[3] [0]; pi[0][i] = ((point[t2][1])*(1.0 — BB[3][i]) —(point[t2+1] [1])eBB[3] [2])/BB[3] [0]; pi[O][2] = ((point[t2][2])*(1.0 — BB[3][1]) —(point[t2+1] [2])*BB[3] [2])/BB[3] [0]; } else{ pi[0][0] = point[t2—i][0]; pi[O][i] = point[t2—1][1]; pi[0] [2] = point [t2—1] [2]; } if(t2 == a — 2){ Appendix M. Codelisting ][O] 3 pi[ = 182 ((point[t2+1][O])*(1.O — bb2) (point[t2] [O])*bbl — (point[t2—1] [O])tbbO)/bb3; — ][i] 3 pi[ = ((point[t2+i][1])s(1.O — (point[t2] [i])*bbi — (point[t2—i] [1])*bbo)/bb3; ][ 3 pi[ ] 2 = — bb2) ((point[t2+1][2])*(1.O — (point[t2] [2])*bbi — (point[t2—1] [2])*bbo)/bb3; — bb2) } else{ pi[3] [0] = point [t2+2] [0]; pi[3] [1] = point [t2+2] [1]; piE3] [2] = point [t2+2) [2]; } } { else if(t2 == 0){ pi[0] [0]((point[t2] [O])*2—(point[t2+1] [0])); piE0] [1]((point[t2] f1])*2—(point[t2+1] [1])); piE0] [2]((point[t2] [2])s2—(point[t2+1] [2])); } { else piEO] [0]point[t2—i] [0]; pi[0] [i]pointft2—i] [1]; pi[O] [2]point[t2—1] [2]; } if(t2 = n — 2){ pi[3] [0]((pointtt2+1] [0])s2—(point[t2] [0])); pi[3] [1]((point [t2+1] [1] )*2—(point [t2] [1])); pi[3] [2]((point Et2+1] [2] )*2—(point [t2] [2])); } else { pi[3] [0]point[t2+2] [0]; pi[3] [1]point[t2+2] [1); pi[3] [2]point[t2+2] [2]; } } p1 [1] [0] point [t2] to]; Appendix M. Codelisting pi[2] [O]point[t2+1] [0]; pi[i] [l]point[t2] [1]; pi[2] [l]=point[t2+l] [1]; pill) [2]9oint[t2] [2]; pi[2) [2]=point[t2+l] [2]; cc = tt*BB; c2 = cc*pi; p[0]c2[O] [0]; p[l]c2[0] [1]; p12] c2 [0] [2]; return p; } vector Curve: :tangent (double t) { tnt t2; double t3 ,bbo ,bbl ,bb2 ,bb3; vector tang; satrix pi(4,3),tt(l,4),c2(l,3),cc(l,4); bbO = BB[0] [O]+BB[l] [0]+BB[2] [0]+BB[3] [0]; bbl = BB[0] [l]+BB[1] [l]+BB[2] [l]+BB[3] [1]; bb2 = EBb] [2]+BB[1] [2]+BB[2] [2]+BB[3] [2]; bbS = EBb] [3]+BB[1] [3]+BB[2] [3]+BB[3] [3]; t2=(int)floor(t); t3t-(double)t2; if(t == n — 1.0; t3 = t2 = it — 2; } tt [0] [0] 3. 0*t3ets; tt[0] [l]2.0*t3; tt[0] [2]l.0; tt[0] [3]0.0; if(aplinetype if(t2 == == pi[O] [0]potnt[t2] [0]; pi[0] [l]point[t2] [1]; pi[0] [2]point[t2] [2]; 183 Appendix M. Codelisting 184 } else{ pit0] [0]point[t2—1] to]; pit0] [1]—point[t2—i] Li]; pi[O] [2]point[t2—i1 [2]; } if(t2 n == — 2){ pit3] [0]pointtt2+i] tO]; p1(3] [i]point[t2+i] [1]; pit3] [2]point[t2+i] [2]; } { else p1(3] [O]point[t2+2] tO]; pi[3] [i]=point[t2+2] Li]; p1(3] [2]point[t2+2] [2]; } } else if (splinetype if(t2 == O){ == pi[O][O] = ((point[t2][O])*(i.O — BB[3][i]) —(point[t2+i] [O])*BB[3] [2])/BB[3] [0]; pi[O][i] = ((point[t2][i])*(i.o — BB[3][i]) —(point [t2+i] [1] )*BB[3] (2]) /BB [3] [0]; pi[O][2] = ((point[t2][2])e(i.0 — BB[3][i]) —(point(t2+i] [2])*BB[3] [2])/BB[3] [0]; } else{ pitO] [0] = p1(0] [1] = point (t2—i] [1]; p1(0] (2] = point [t2—i] [2]; point (t2—i] (0]; } if(t2 = n — pi[3][O) 2) = { ((point[t2i-1]tO])s(i.O — (point [t2] [0] )*bbi — (point(t2—i] [0])sbbO)/bb3; pi[a](i] = ((point[t2+i](i])s(i.0 — (point[t21 [i])*bbi — (point[t2—i] (i])*bbo)/bb3; pi(3](2] = ((point(t2+i](2])s(i.0 — bb2) — bb2) — bb2) Appendix M. Codelisting 185 — (point[t2] [2])*bbl — (point[t2—l][2])*bbo)/bbs; } else{ p113] [0] = point [t2+2] 10]; p113] [1] = point [t2+2] [1]; p113] [2] = point [t2+2] 12]; } } else{ if(t2 == 0){ p1 [0] [0] ( (point [t2] [0]) *2—(point [t2+l] [0])); p1 [0] [1] ( (point [t2] [1] )s2—(point [t2+l] [1])); p110] [2] ( (point [t2] [2] )*2—(polnt [t2+l] [2])); } else { p110] [0]point [t2—1] [0]; p110] [l]point[t2—l] It]; p110] [2]point[t2—l] [2]; } if(t2 == n — 2){ p113] [0]((point[t2+t] [0])*2—(point[t2] [0])); pi[3] [t]((point[t2+l] [1])*2—(point[t2] [1])); pi [3] [2] ( (point 1t2+l] [2]) *2— (point [t2] [2])); } else { p113] [0]point[t2+2] [0]; p113] [l]point[t2+2] [1]; p113] 12]point[t2+2] [2]; } } pill] [0]point[t2] [0]; p112] l0]pointlt2+l] [0]; pill] ll]point[t2] Ii]; p112] ll]pointlt2+l] [1]; pill] 12]pointlt2] [2]; p1(2] 12]point[t2+1] [2]; cc = tt*BB; c2 = cc*pi; Appendix M. Codelisting tang[0)c2[0] [0]; tang [1] c2 [0] [1]; tang [2] c2 [0] [2]; return tang; } vector Curve: :curvature(double t) { mt t2; double t3,bb0,bbl ,bb2,bb3; vector curv; natrix pi(4,3),tt(1,4),c2(1,3),cc(1,4); bb0 = BB[0] [0]+BB[1] [0]+BB[2) [0]+BB[3] [0]; bbl = BB[0] [1]+BB[1] [1]+BB[2] [1]+BB[3] [1]; bb2 • BB[0] [2]+BB[1] [2]+BB[2] [2]+BB[3] [2]; bb3 = BB[0] [3]+BB[1] [3]+BB[2] [3]+BB[3] [3]; t2(int)floor(t); t3t—(double)t2; i±(t t3 = 1.0; t2 = n — == n — 2; } tt[0] [0]6.0st3; tt[O] [1]2.0; tt[O] [2]0.0; tt[0] [3]0.0; if(splinetype == ‘1’) { if(t2 == o){ pi[O] [0]point[t2] [0]; pi[0] [1]point[t2] [1]; pi[0] [2]point[t2] [2]; } else { pi[0] [0]point[t2—1] [0]; pi[0] [1]point[t2—1] [1]; pi[0] [2]point[t2—1] [2]; } 186 Appendix M. Codelisting if(t2 == n — 187 2){ pi[3] [OF’point[t2+1] [0]; pi[3] [1]point[t2+i] [1]; pi[3] [2]point[t2+i] [2]; } { else pi[3] [Ofrpoint[t2+2] [0]; pi[a] [1]point[t2+2] [1]; pi[3] [2]point[t2+2] [2]; } } else if (splinetype == ‘2’) { if(t2 == 0) { pi[0J[0] = ((point[t2][0])*(i.0 — BB[3][i]) —(point [t2+1] [0])sBB[3] [2])/BB[3] [0]; pi[0][t] = ((point[t21[1])*(1.0 — BB[3][1]) —(point[t2+i] [1])sBB[3] [2))/BB[3] [0); pi[0][2] = ((point[t2][2])*(1.0 — BB[3][1]) —(point [t2+1] [2])CBB[3] [2])/BB[3] [0]; } else { pi[0][0] = point[t2—1][0]; pifO] [1] = point [t2—i] [ii; pi[0][2] = point[t2—1][2]; } if(t2 == n — 2) { pi[3)[O] = ((point[t2+1][0])*(i.0 — (point [t2) [0))sbbi — (point[t2—1] [0])ebbo)/bb3; pi[3][1] = ((point[t2+i][i])*(1.0 — (point[t2] [1])sbbi — (point[t2—1] [1))sbbO)/bba; pi[3][2] — = ((point[t2+i)[2))*(1.0 (point[t2] [2])*bbi — bb2) — bb2) — bb2) Appendix M. Codelisting — (point [t2—1) [2) )*bbO)/bb3; } else { pi[3] [0] = point [t2+2) [0]; pi[3] [1] = point [t2+2] [1]; pi[3][2] = point[t2+2][2]; } } else { if(t2 == 0){ ( (point [t2] [0]) *2— (point [t2+1] [0])); pi[0] [0] pi[0] [1] = ((point [t2] Li] )*2—(point [t2+i] [1))); pi[0] 12]((pointEt2] [2))*2—(point[t2+i] [2])); } { else pi[O] [0)point[t2—i] [0]; pi[0] [i]point[t2—i] [ii; pi[0] [2]point[t2—i] [2]; } == if(t2 n — 2){ pi[3] [0)((point[t2+i] [0])*2—(point[t2] [0])); pi[3] [i]( (point [t2+1) [1)) *2— (point [t2) Li])); pi[3) [2] = ((point [t2+1] [2)) *2— (point [t2] [2))); else { pi[3] [0]point[t2+2] [0]; pi[3] [i]point[t2+2] [1); pi[3] [2]point[t2+2] [2); } } pi[i] [0]point [t2] [0]; pi[2] [0]point[t2+i] [0); piLl] [l]point[t2] [1]; 188 Appendix M. Codelisting 189 pi[2) [1)i:point[t2+1) [1]; piLl) [2frpoint [t2) [2); pi[2) [2)=point[t2+1) [2); cc = tt*BB; c2 = cc*pi; cnn 10) c2 [0) [0); cun[1)=c2[0) [1); cun[2)c2[0) [2); return curT; } double Curve: :deriv_2c(double t, mt 1) { jut t2; double t3,deriv; natrix tt(l,4),cc(1,4); t2(int)floor(t); t3t- (double)t2; tf(t t3 = 1.0; t2 = a — == a — 2; } if( i C (t2—1) II i > (t2+2) ) return 0.0; tt[0) [0]3.0*t3.t3; tt[0) [1]2.0et3; tt[0] [2fr1.0; tt[0] [3fr0.0; cc = tt*BB; if( splinetype == ‘2’) { if(t == t2 && t 0.0 U t ! n—2) { cc[0) [0]cc[0) [1); cc[0) [2]cc[0) [3); } if(i == t2 — 1) deny = cc[0][0]; Appendix M. Codelisting else 11(1 t2) den, == = 190 cc[O][1]+(t20.O?2.0*cc[O][O]:O.0)+ (t2n—2?—cc[0] [3] :0.0); I/Was t2n—2, Bug caught Aug 8,1992 else if(i t2 + i) deny = cc[0][2]+(t20.0’?—cc[0][0]:0.0)+ (t2n—2?2.O*cc[0] [3] :0.0); I/Was t2n—2, Bug caught Aug 8, 1992 else if(i t2 == else deny + 2) deny = cc[0][3]; 0.0; = } else if( splinetype == ‘1’) { jf( 1 t2 = — 1) deny = cc[0][0]; = cc[O][1] +(t20.O?cc[O][0]:O.O); else if(i = t2) deny else if(i == t2 + 1) deny = cc[0][2] itCi == t2 + 2) deny = cc[0][3]; else else deny + (t2n—2?cc[0] [3] :0.0); 0.0; = } else { if(t == t2 U t 0.0 U t n-2) { cc[0] [0]cc[0] [1]; cc[0] [2]cc[0] [3]; } i:f( ± 12 == else if(i — == 1) deny = cc[0][0]; t2) deny = cc[0][1]+(t2flO.0?2.Oscc[0][0]:0.0) +(t2n—27—cc[0] [3] :0.0); else if(i == t2 + 1) deny = cc[0][2]+ (t20.0?—cc[0] [0] :0.0) +(t2n—2?2.O*cc[0] [3] :0.0); else ±1(1 else deny == = t2 + 2) deny = cc[0][3]; 0.0; } return deny; } double Cunve::deniv..c(double t, mt 1) { jut t2; Appendix M. Codelisüng 191 double t3,deriv; atrix tt(1,4),cc(1,4); t2(int)floor(t); t3t—(double)t2; if(t t3 = t2 = == n — 1.0; — 2; } if( i < (t2—1) II i > (t2+2) ) return 0.0; tt [0] [0] t3*t3*t3; tt[0] [1]t3*t3; tt[0] [2]t3; tt[0] [3]1.0; cc = tt*BB; if( splinetype == ‘2’) { if(i t2 1) deny — cc[0)[0); = else if(i t2) deny else if(i t2 + 1) deny = cc[0J[2] t2 + 2) deny = cc[0][3); else if(i == else deny = cc[0][1] +(t20.0?2.0*cc[O) [0] :0.O)+(t2u—2?—cc[O] [3] :0.0); + (t20.0?—cc[0)[0]:0.0)+(t2u—2?2.0*cc[0][3]:Q.O); 0.0; = } else if( splinetype == ‘1’) { if( i == t2 — else if(i 1) deny = cc[0][0J; t2) deny = cc[0][1] +(t2O.O?cc[Q][O]:O.O); else if(i == t2 + 1) deny = cc[0J[2] else if(i == t2 + 2) deny = cc[0][3); else deny + (t2n—2?cc[0] [3] :0.0); 0.0; = } else { if( i t2 — 1) den, = cc[0][0]; = cc[0][1] +(t2O.O?2.O*cc[O][O]:O.O)+(t2—2?—cc[O][3]:O.O); else if(i == t2) den, else if(i == t2 + 1) deny = cc[0J[21 else if(i == t2 + 2) deny = cc[0][3]; else deny } = 0.0; + (t20.0?—cc[0]f0):0.0)+(t2u—2?2.0*cc[0][3]:0.0); Appendix M. Codelisting return dent; } void Curve: :initbasis(char spitype, double bi, double b2) { double del; if(spltype ‘1’) { -bi; BB[0][0) 2.0—bi; BB[0][1) BB[0][2] bl—2.0; = bi; BB[0][3] BB[1] [0] = 2.0*bl; BB[1][1] = bl—3.0; BB[1][2] = 3.0—2.0*bl; BB[1][31 = —bi; —bi; BB[2][0] BB[2][1] = 0.0; BB[2][2] = bi; BB[2][3] 0.0; BB[3][0] 0.0; BB[3][1] = 1.0; BB[3][2] = 0.0; BB[3][3J = 0.0; } else if(spltype == ‘2’) { del = (b2)+(2.0*(bl)s(bl))+(4.0*(bl)*(bl))+(4.0*(bl))+2.0; BB[0] [0] = —2.0*(bl)*(bl)*(bl)/del; BB[0][1) = 2.0*((b2)+((bl)*(bl)*(bl))+((bl)*(bl))+(bl))/del; BB[0] [2] = —2.0*((b2)+(bl)*(bl)+(bl)+1.0)/del; BB[0][3] 2.0/del; BB[1][0] = BB[1][1] = 6.0*((bl)*(bl)*(bl))/del; 3.0*((b2)+2.0*((bl)*(bl)*(bl))+2.0*((bl)*(bl)))/del; 3.0*((b2)+2.0*((bl)*(bl)))/del; BB[1][2] BB[1][3] = 0.0; BB[2][0] = —6.0$((bl)*(bl)*(bl))/del; 192 Appendix M. Codelisting BB[2][1] = 6.0*((bl)*(bl)*(bl)—(bl))/del; BB[2][2] = 6.05(M)/del; BB[2][3] = 0.0; BBI3][O] = 2.0S(bl)*(bl)*(bl)/del; BB[3]11] = ((b2)+4.0s((bl)s(bl))+4.0*(bl))/del; BB[3]12] = 2.0/del; BB[3][3] = 0.0; } else { del = (b2)+(2 .0*(bl)*(bl))+(4.0*(bl)*(bl))+(4.0s(bl))+2 .0; BB[0] 10) = Ii) = 2.0t((b2)+((bl)*(bl)*(bl))+((bl)*(bl))+(bl))/del; BB[0] 12] = —2.0*((b2)+(bl)*(bl)+(bl)+1.0)/del; EB[0] 2.0*(bl)*(bl)*(bl)/del; BBEO][3] = 2.0/del; BB11]10] = 6.0*((bl)*(bl)*(bl))/del; BB11][1] = BBI1)[2) = 3.0*((b2)+2.0*((bl)S(bl)))/del; BB[1][3] = 0.0; 3.0*((b2)+2.0*((bl)*(bl)*(bl))+2.0*((bl)*(bl)))/del; BB[2][0] = —6.0S((bl)*(bl)*(bl))/del; BB[2][1] = 6.0*((bl)*(bl)*(bl)—(bl))/del; BB[2][2] = 6.0*(bl)/del; BB[2][3] = 0.0; BB[3][0] = 2.0*(bl)s(bl)*(bl)/del; BB[3][1] = ((b2)+4.0*((bl)*(bl))+4.0*(bl))/del; BEla] 12] = 2.0/del; BB[3][3] = 0.0; } BBBBB; } ofstreea& Filelisto: :operatorlj (hit a) { Filelisto ste.pthis; for(; a>0; a-—) teap return te.p—>file; } = te.p—>aext; 193 Appendix M. Codelisting ifstreaat Filelisti: :operatorfl (tnt 194 11) { Filelist i ttespthis; for(; n>O; n--) te = te.p->next; return teap—>file; } M.1.9 ode.cpp tiuclude “ode.h” inline double abs(double a){ return a<O.O?-a:a; } dynamic_systea: :dynaaic_systea(vector (sf)(double, vectort), double start_time, vectort initial_state, vector terrors) state (initial_state) ,error_scale (errorsO?init ial_state: terrors) { t iaentart _t iae; step_sizeO.O; if (errorsO){ for (tat isize(state)—1;i>0;i——) error_scale[i] = 1.Oe-4; } derivativef; } doublet dynamic_systea: : when() { return time; } void dynaaic_systea: :reset() { tiaeO.O; } doublet dyneaic_systea: :operatorfl(int i) Appendix M. Codelisting { return state[i]; } vector dynamic_system: :operatorO(double new_time) { if(step_sizeo .0) step_sizeabs(new_time—tiae)/2 .0; rk(time,step_size,new_tiae,state,derivative,error_scale); timenew_tiae; return(state); } vector dynamic_system: :step(double delta) { if(step_sizeo .0) step_stzedelta/2 .0; rk(time,step_size,time+delta,state,derivative,error_scale); time + delta; return state; } 195 Index Aerospace, 2 Computer Language Chosen, 83 Alignment of End Generators, 33 concept of mobility, 38 alignment process, 36 Conical Type Surface, 77 Allowment of Different Number of Con Constraints Defining Modified Conven trol Vertices, 48 tional Approach (Modified Nolan’s Arctic Fishing Vessel, 100 Approach) in Order to Create a Areas of Application, 1 Normal Directrix AutoCAD, 96 lan’s Approach) in Order to Create a Normal Directrix, 42 Barsky, 17 continuity, 18 Basis functions, 16 Controlling Alignment with Mobility, 39 C++ Classes, 85 Creation of a Normal Directrix from Two CAD Program Environment, 96 Space Curves, 30 CAD Program Selected, 96 Demonstration Examples, 98 Change in Angle With Respect to Unit dependent parametric variables, 73 Motion of a Control Vertex, 34 Derived Equations Yielding Intersections Class Curve, 91 of Surfaces, 74 Class matrix, 90 develop.cpp, 179 Class ODE, 93 develop.h, 162 Class Surface, 91 Developable Mobius Strip, 98 Class vector, 85 developable mobius strip, 39 Clement, 9 developable surface, 4 Code Portability to Different Platforms, directrix, 4 85 196 Index 197 Downhill Simplex Method, 67 matrix.h, 161 Dunwoody, 10 Methodologies, 5 Dunwoody and Konesky, 31 Mobius Strip, 38 Equations Yielding Approximated Nor Modern Approach Utilizing a Single Nor mal Directrix, 31 mal Directrix, 52 explicit non parametric form, 14 F117A stealth fighter, 2 fanning, 49 flat plate equation, 81 Flat Plate Layout, 80 flat plate layout, 81 Naval Architecture, 1 new approach, 12 Nolan, 8 non parametric implicit formulation, 15 non-equal number of control vertices, 46 non-uniform curves, 18 Normal Directrix Approach, 10 geometric continuity, 19 normal directrix control vertices, 50 in-plane curvature, 46 ode.cpp, 194 independent parametric variable, 73 ode.h, 164 Integral Chosen to Minimize, 64 Offset of Normal Directrix, 47 Interplate Angles, 80 OOA Intersection of Developable Surfaces, 73 OOD Intersection of Developable Surfaces and OOP Object-Oriented Programming, 84 Flat Plate Layout, 73 Kilgore’s Manual Graphical Solution, 6 - - Object-Oriented Analysis, 84 Object-Oriented Design, 84 - Open Architecture, 96 Optimization Function, 60 Optimization of a Normal Directrix, 60 least squares approximation, 64 out-of-plane curvature, 46 Manufacturing, 2 parametric continuity, 18 matrix.cpp, 169 parametric curves, 15 Index 198 PC 386/486 Environment, 95 Uniform Non-Rational Beta-Spline, 24 phantom, 78 Uniform Non-Rational Tension Catmull Phantom Surfaces, 78 Rom Spline, 21 phantom surfaces, 77 Utilizing Modified Conventional Approach PNA, 6 to Approximate a Single Direc rate-of-rotation of a generator with re spect to motion of the control ver tices along the normal directrix ntrol vertices along the normal directrix, 34 rational continuity, 20 ruled surface, 3 - Simple Conical Developable, 100 Simplex Parameters Used, 64 small in and out-of-plane curvature, 41 solve.cpp, 175 Splines, 14 Textiles, 3 The Downhill Simplex Method, 68 threshold, 42 tolerance, 42 Types of Splines Chosen, 14 UBC Series Fishing Vessel, 100 uniform curves, 18 trix Directrix, 50 vector.cpp, 165 vector.h, 160
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Computer aided design of developable surfaces
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Computer aided design of developable surfaces Konesky, Brian E. 1993
pdf
Page Metadata
Item Metadata
Title | Computer aided design of developable surfaces |
Creator |
Konesky, Brian E. |
Date | 1993 |
Date Issued | 2009-02-23T20:02:34Z |
Description | The design of objects employing developable surfaces is of engineering importance because of the relative ease with which developable surfaces can be manufactured. The problem of designing developable surfaces is not new. Two space curves, defining the edges of the surface, are first created, then a set of rulings are constructed between the space curves under the constraint of developability. A problem with existing algorithms for designing developable surfaces is the tendency to introduce non developable portions of the surface; areas of regression. A more reliable solution to the problem of creating a developable surface is proposed. The key to the method is to define the developable surface in terms of a normal direc trix. The shape of the normal directrix defines the shape of the developable surface. Algorithms are defined to compute the shape of a normal directrix from a pair of space curves. A non-linear optimization technique was implemented to further refine the shape of the developable surface, but failed to yield satisfactory results. Other algorithms were also created to intersect adjacent developable surfaces and to generate the flat plate lay outs. The algorithms were implemented using the C++ programming language and the AutoCAD CAD package. Recommendations for further work are given. The design of objects employing developable surfaces is of engineering importance because of the relative ease with which developable surfaces can be manufactured. The problem of designing developable surfaces is not new. Two space curves, defining the edges of the surface, are first created, then a set of rulings are constructed between the space curves under the constraint of developability. A problem with existing algorithms for designing developable surfaces is the tendency to introduce non developable portions of the surface; areas of regression. A more reliable solution to the problem of creating a developable surface is proposed. The key to the method is to define the developable surface in terms of a normal directrix. The shape of the normal directrix defines the shape of the developable surface. Algorithms are defined to compute the shape of a normal directrix from a pair of space curves. A non-linear optimization technique was implemented to further refine the shape of the developable surface, but failed to yield satisfactory results. Other algorithms were also created to intersect adjacent developable surfaces and to generate the flat plate lay outs. The algorithms were implemented using the C++ programming language and the AutoCAD CAD package. Recommendations for further work are given. |
Extent | 3735080 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-02-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080847 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1994-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/4930 |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- ubc_1994-0078.pdf [ 3.56MB ]
- [if-you-see-this-DO-NOT-CLICK]
- Metadata
- JSON: 1.0080847.json
- JSON-LD: 1.0080847+ld.json
- RDF/XML (Pretty): 1.0080847.xml
- RDF/JSON: 1.0080847+rdf.json
- Turtle: 1.0080847+rdf-turtle.txt
- N-Triples: 1.0080847+rdf-ntriples.txt
- Original Record: 1.0080847 +original-record.json
- Full Text
- 1.0080847.txt
- Citation
- 1.0080847.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
Russia | 25 | 0 |
China | 23 | 2 |
Canada | 9 | 4 |
Germany | 9 | 2 |
France | 8 | 0 |
United States | 8 | 4 |
Romania | 1 | 0 |
United Kingdom | 1 | 0 |
Unknown | 1 | 0 |
Indonesia | 1 | 0 |
Ukraine | 1 | 0 |
Turkey | 1 | 0 |
Poland | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 40 | 1 |
Beijing | 22 | 0 |
Wadgassen | 8 | 0 |
Ashburn | 4 | 0 |
Vancouver | 4 | 0 |
Hamburg | 1 | 0 |
Denver | 1 | 3 |
Shenzhen | 1 | 2 |
Ankara | 1 | 0 |
Suceava | 1 | 0 |
Overland Park | 1 | 0 |
Montreal | 1 | 2 |
Buffalo | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080847/manifest