COMPUTER AIDED DESIGN OF DEVELOPABLE SURFACESByBrian E. KoneskyB.A.Sc. (Mechanical Engineering) University of British ColumbiaA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE STUDIESMECHANICAL ENGINEERINGWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1993© Brian E. Konesky, 1993In presenting this thesis in partial fulfilment of the requirements for an advanced degree atthe University of British Columbia, I agree that the Library shall make it freely availablefor reference and study. I further agree that permission for extensive copying of thisthesis for scholarly purposes may be granted by the head of my department or by hisor her representatives. It is understood that copying or publication of this thesis forfinancial gain shall not be allowed without my written permission.Mechanical EngineeringThe University of British Columbia2075 Wesbrook PlaceVancouver, CanadaV6T 1W5Date:AbstractThe design of objects employing developable surfaces is of engineering importance becauseof the relative ease with which developable surfaces can be manufactured. The problemof designing developable surfaces is not new. Two space curves, defining the edges of thesurface, are first created, then a set of rulings are constructed between the space curvesunder the constraint of developability. A problem with existing algorithms for designingdevelopable 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 spacecurves. A non-linear optimization technique was implemented to further refine the shapeof the developable surface, but failed to yield satisfactory results. Other algorithms werealso created to intersect adjacent developable surfaces and to generate the flat plate layouts. The algorithms were implemented using the C++ programming language and theAutoCAD CAD package. Recommendations for further work are given.11Table of ContentsAbstractList of Tables ixList of Figures xNomenclature xiiiAcknowledgement xvii1 Introduction1.1 Areas of ApplicationNaval ArchitectureAerospaceManufacturingTextiles1.2 Definitions and Terminology1.3 Methodologies1.3.1 Kilgore’s Solution1.3.2 Nolan and Clement Computer Approach .1.3.3 Normal Directrix Approach by Dunwoody .1.4 New Approach Research Objectives2 Splines 14111223355810121112.1 Types of Splines Chosen.142.2 Uniform Parametric, Geometric and Non-Rational Continuity Requirements 182.3 Uniform Non-Rational Tension Catmull-Rom Spline 212.4 Uniform Non-Rational Beta-Spline 243 Creation of a Normal Directrix from Two Space Curves 303.1 Modern Approach Utilizing a Single Normal Directrix 313.1.1 Constraints Governing Modern Approach 323.1.2 Alignment of End Generators 33Change in Angle With Respect to Unit Motion of a Control Vertex 34The Alignment Process and the Concept of Mobility 363.1.3 Analysis of Results 38Mobius Strip Demonstrating Flexibility 38Controlling Alignment with Mobility 39Problems arising when Surface has Small In and Out-of-plane Curvature 413.2 Constraints Defining Modified Conventional Approach (Modified Nolan’sApproach) in Order to Create a Normal Directrix 423.2.1 Offset of Normal Directrix 473.2.2 Allowment of Different Number of Control Vertices 483.2.3 Present Problems with Current Solution and Improvisation Implemented 493.3 Utilizing Modified Conventional Approach to Approximate a Single Directrix 503.3.1 Equations Yielding Approximated Normal Directrix 523.3.2 Results and Present Problems 53iv4 Optimization of a Normal Directrix 604.1 Objective 604.2 Optimization Function 604.2.1 Integral Chosen to Minimize and Simplex Parameters Used . 644.2.2 Present Problems 664.3 Downhill Simplex Method 674.3.1 Explanation of The Downhill Simplex Method . . 684.3.2 Results and Present Problems 694.4 Examples 694.4.1 Example of Single Surface with Conical Properties 694.4.2 Example of UBC Series Demonstrating Present Problems 705 Intersection of Developable Surfaces and Flat Plate Layout 735.1 Intersection of Developable Surfaces 735.1.1 Derived Equations Yielding Intersections of Surfaces 745.1.2 Results and Present Problems 77Conical Type Surface Tested 77Criterion of Phantom Surfaces 785.2 Flat Plate Layout 805.2.1 Derived Equations giving Flat Plates Dimensions and InterplateAngles 805.2.2 Format of Output 816 Choice of Computer Language and CAD Program 8361. Computer Language Chosen 836.1.1 OOP - Object Oriented Programming 846.1.2 Portability to Different Platforms 85V6.1.3 C++: Classes and FeaturesClass vectorClass matrixClass Curve and SurfaceClass ODE6.2 Computer Platform Selected6.2.1 PC 386/486 Environment6.2.2 Computer Language Selection in this Platform6.3 CAD Program Selected6.3.1 CAD Program Environment and Open Architecture6.3.2 Present Limitations85859091939495959696977 Demonstration Examples7.1 Developable Mobius Strip7.2 Simple Conical Developable7.3 Arctic Fishing Vessel7.4 UBC Series Fishing Vessel8 Conclusions and Recommendations8.1 ConclusionsNormal Directrix from Two Space CurvesIntersection of Developable Surfaces and Flat Plate Layout . .Implementation8.2 RecommendationsBibliography 108Appendices 1109898100100100105105105106106107viA Mathematical Notation for Partial Differentiation 110B Derivation of Developable Surface 111B.1 Constraints which Define the Developable Surface 111B.2 Proof Using Constraints 113C Derivation of Rate of Rotation of Generator Differential Equation 115D Tension Catmull-Rom Spline 122D.l Phantom point, P(O), at 0 123D.2 Phantom point, P(1), at n 125E Beta-Spline 127E.l Phantom point, P(0), at 0 128E.2 Phantompoint,P(1),atn 131F Derivation of Normal Directrix Control Vertices 134G Modified Conventional Approach Derivation 139H Relating the Space Curves to the Developable Surface 142I Intersection of Developable Surfaces 146J Derivation of Flat Plates 151K O.D.E. Class, Adaptive Step Size and Runge-Kutta Method 153L Non-Linear Optimization Downhill Simplex Method 155M Codelisting 160vi’M.1 Class Tools Used.160M.1.1 vector.h 160M.1.2 matrix.h 161M.1.3 develop.h 162M.1.4 ode.h 164M.L5 vector.cpp 165M.1.6 ma.trix.cpp 169M.L7 solve.cpp 175M.1.8 develop.cpp 179M.1.9 ode.cpp 194viiiList of Tables5.1 Flat plate data.82L.1 Values of variables used for Simplex Method 157ixList of Figures1.1 Developable Surface Single Chine Hull 21.2 Manual method of generating developable surfaces 31.3 F117A stealth fighter 41.4 Definitions of Terminology of a Developable Surface 51.5 Kilgore’s Method of Creating a Developable Surface 61.6 Kilgore’s Manual Graphical Solution 71.7 Nolan’s Vector Description of Computer Solution 91.8 Normal Directrix Approach by Dunwoody 112.1 First Three Levels of Parametric Continuity 192.2 Comparing Geometric and Parametric Continuity 212.3 Catmull-Rom Splines 252.4 Beta Splines 293.1 Derivation of Developable Surface 333.2 Vector Locations and Corresponding Angles 343.3 Control Vertices and End Phantom Point 373.4 Robustness of Modern Approach Showing Mobius Strip 393.5 Mobius Strip 403.6 Modern Approach Showing Full Mobility of All Points 413.7 Provision Made for When Surface Very Near Flat 423.8 Relating Modified Conventional Approach Space Curves and Normal Directrix 43x3.93.103.113.123.133.143.153.163.173.184.14.24.34.44.64.7Developable Surface FailureDevelopable Surface success, tolerance 0.001 .44485053545556575859616567687071725.15.25.35.45.574• . . . 78• • . • 7980816.1 Computer Platform Selected Initially7.1 Developable mobius Strip9499In-plane and out-of-plane componentsPositional Offset of Normal DirectrixFanning Occurring When Developable Surface Nearly FlatConic-Like Section with Flatness Tolerance Specified at 0.001Conic-Like Section With No Flatness Tolerance SpecifiedDevelopable Surface Procedure Comparison, tolerance 0.00 1Developable Surface success, tolerance 0.01Developable Surface Procedure Comparison, tolerance 0.01Developable Surface success, tolerance 0.1Developable Surface Procedure Comparison, tolerance 0.1Relating the Distance Between Space Curves and Surface . .Area calculated under space curvesPossible inherent failure due to certain geometryAnalogy of a Simplex For The Downhill Simplex Method . .4.5 Conical-like developable surfaceIntersection of Three Developable SurfacesConical-type Surfaces IntersectionsConical-type Surfaces IntersectionsFlat plate layout derivationDevelopable surface and flat plate layoutxi7.2 Conical-type Surfaces Intersections 1017.3 Arctic Vessel Conventional Approach 1027.4 Arctic Vessel Modern Approach 1037.5 UBC series Vessel Conventional Approach 104B.1 Derivation of Developable Surface 111B.2 Derivation showing normal is invariant along a generator . . 112C.l Vector locations and corresponding angles 115G.l Orientation of Space Curves and Directrix 1391.1 Intersection of Three Developable Surfaces 146L.l Analogy of a Simplex for the Downhill Simplex Method 156xiiNomenclature45j,j 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 Describing Directrix.Rate of Change of Angle of Generator With Respect to Motion of Parameter Describing Directrix.In Plane Rotation of G’ With Respect to G.a Parameter Describing Distance Along Generator.4a 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.XIIICurvature Vector Along Right Adjacent Space Curve.. The Developable Surface Generator Vector.dg:ri Change of Position of Generator With Respect to Motion of Parameter DescribingDirectrix.d2genRate of Change of Generator With Respect to Motion of Parameter DescribingdadtDirectrix.dgenThe 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 Generatorxivq(t, s) Point on Surface Determined by Parametric Values t and Scalar st Independent Parametric Variable Describing Position Along Directrix.T Unit Directrix Tangent Vector.Change of Position of Unit Tangent With Respect To Motion of Parameter Describing 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.xvA. EinsteinKeep it simple:as simple as possible,but no simplerENGINEERINGThe scientist analyzes what is.The engineer creates what has never been.The engineer scientist analyzes what isImagines what should beCreates what has never beenAnalyzes the results of the creation.By Gunnar SchieniusxviAcknowledgementI would like to express my sincere thanks to Dr. A. B. Dunwoody for his enthusiasticsupport, supervision and direction on this project.I would also like to acknowledge the financial support from the NSERC of Canadawhose financial support made this project possible.xviiChapter 1Introduction1.1 Areas of ApplicationDevelopable surfaces form a special class of surfaces which are very useful in many practical situations. Developable surfaces have many applications. A few applications arecited below.Naval Architecture Up until recently in history ship hulls were made of various typesof 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 structureof small pieces of wood. If a skin of sheet material can be designed for low labour cost inconstruction, simple tools, and economy in repair, its engineering superiority and eventualeconomic advantage make it at once preferable to planks” [15]. To lower the labour costin 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 withthat of the best possible in compound hull construction [15]. Making the constructionsurface 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 surfacehulls for fishing vessels.1Chapter 1. Introduction 2Aerospace Today’s modern aircraft are more complex in design and cost savings inmanufacturing is crucial. Aircraft fuselages, wings and other smaller components canbe produced by developable surfaces if considered in the design stage. One of the mostrecent “high-tech” declassified aircraft is the “wobblin-gobblin”, ie. the “Fl 1 7A”, (seefigure 1.3) which has a low radar profile signature by using flat plates, is yet anotherexample of using developable surfaces in an ingenious manner.Manufacturing In the area of manufacturing two new areas involving the developability criterion are in the application of peripheral milling and rolling. Each pass of anend mill cutting peripherally follows a developable surface.Figure 1.1: Developable Surface Single Chine HullChapter 1. Introduction 3Figure 1.2: Manual method of generating developable surfacesAnother area is in the application of rolling. Developability must also hold for thisprocess.Both of these applications will not be discussed but only cited as other exampleswhere developable surface criterion must hold due to the geometry of the application.Textiles Another example of developability is in the textile industry. Clothing ismanufactured 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 TerminologyThe developable surface is a subset of the class of ruled surfaces. The definitions of aruled 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 4Figure 1.3: F117A stealth fighterwhose direction is determined by successive values of a parameter, moving continuouslyalong 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 surfacehaving the same tangential plane on one and the same generator” [15].Figure 1.4 below illustrates the terms introduced similar to the definitions presentedby G.D. Aguilar[2j. The directrix must lie in the surface and that each plane tangentialto the surface must also be tangential to the directrix.It must also be noted that with two space curves a ruled or compound surface mayexist but no developable surface may be possible. Another theorem should also be notedthat, “If two space curves lie in any developable surface, they lie in one and only one suchsurface” [15]. If the generators do not intersect anywhere, then the surface is developable.Chapter 1. Introduction 5Space CurvesFigure 1.4: Definitions of Terminology of a Developable SurfaceIn figure 1.4 the areas where the generators overlap are known as areas of regression. Theboundary of an area of regression outlined in figure 1.4 is called the edge of regression[71.1.3 MethodologiesAll existing methodologies for the design of developable surfaces start with the definitionof two edges of the surface. Then, a set of generators is fit between those edges to definea developable surface.1.3.1 Kilgore’s SolutionThe general method of matching developable surfaces to desired curves is an arduoustask. The method assumes that the developable surface is either conical, cylindrical, ora combination of both. So, one tries a succession of surfaces until one is found to fitGenerator or RulingEdge of RegressionChapter 1. Introduction 6approximately. If the designer cannot find an exact solution his usual solution is to alterthe 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 SurfaceKilgore examined this unique art and proposed a method for direct generation ofdevelopable surfaces from given beginnings. This method provides a manual graphicalsolution of surfaces to fit the curves, rather than to require alteration of the curves to fitthe surfaces.Kilgore’s manual graphical solution is described in his paper[15] as well as a comprehensive 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 fromthe Principles of Naval Architecture[ll]. See Figure 1.6.One can see that manual graphical solutions are only as accurate as the skills andprecision of the naval architect. This poses problems in error of the final solution as toits accuracy.z n -v I ni c-fl 0 ‘1 z > r n x -I rn n -I C rnpçn r p‘ a r-p I II noq I-, CD -a cm 0 WI CD C)‘1 C n 0 S 011 C .4 I, 0 3 0 I 4I-4Chapter 1. Introduction 8This led to implementing a mathematical solution once computers had evolved tothe point where it was a viable alternative. The first computer solutions which had amoderate success were those by Nolan [17] and Clements [7].1.3.2 Nolan and Clement Computer ApproachOne of the first well known published works involving a computer-aided approach todevelopable surface design was by T. J. Nolan [17]. In his paper he emphasized that amathematical approach utilizing a computer proved to yield a substanitial increase inspeed and precision for calculating a developable surface. Nolan noted that an infinitenumber of surfaces can be found to span the curves but the developable surface is uniquein that it requires the minimum strain energy of flexture and that in a developable surfacebending is restricted to nonintersecting axes lying in the surface so that section moduliand bending stresses are minimized for any given radius of curvature. As a result, adevelopable surface can be formed elastically from a plane sheet, while the surface fittingthe same pair of curves but having compound curvature must undergo a costly plasticforming process. Nolan defines a developable surface as, “a developable surface spanninga pair of curves in space may be defined as the locus of straight lines or “rulings” whichrepresent the line of contact of a plane which is tangent to the surface. The rulings areneutral axes of bending and must not intersect within the surface” [17].Nolan’s approach utilizes a Theilheimer spline interpolation and a vector representation to create a mathematical description for a developable surface. The approach isrelatively simple, involving a representation of the tangents, normals and rulings in vector form. He iteratively solves his mathematical model to yield a zero angle for the crossproduct of the two normals of the space curves. See figure 1.7.Nolan’s vector approach calculates the normal of each space curve as N1 = 1? x T1,and N2 = R x T, where R is the ruling and T1 is the tangent calculated at that pointChapter 1. Introduction 9Cuef(uxp(t)eg(v)NiFigure 1.7: Nolan’s Vector Description of Computer Solutionalong space curve 1 and T2 is the tangent calculated at a point on space curve 2. He thenuses one of the constraints of a developable surface, namely that N1 x N2 = 0.This approach is moderately successful in that for very simple surfaces it may yielda resulting developable surface. However, all too often the rulings either cross or “fan”yielding an unrealistic or nonusable surface. Also, the Theilheimer spline interpolation isnot of parametric form resulting in singularities which may occur in a three dimensionalrepresentation of the surface.Clement’s Solution to the problem was published approximately ten years later utilizing cubic spline functions, again in non-parametric form[7]. He states, “Between eachpair of chine lines that ruled surface generated which has the same tangent plane at alldfdunerator9=(;f)TargentPlaneChapter 1. Introduction 10points of each generator or ruling line. A procedure based on the multiconic developmentof a surface is used to modify the given chine lines to ensure that no ruling lines intersect 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 hadproblems arising at singularities and generators crossing in areas of regression on thesurfaces.1.3.3 Normal Directrix Approach by DunwoodyThe approach taken by Dunwoody is unique in that it defines a developable surface bymeans of a normal directrix and an initial generator. The directrix is modelled by aparametric cubic spline; initially a uniform non-rational B-spline was used. Informationabout 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 developable surface. Refering to figure 1.8 to show the geometry, the following constraints areshown and the resulting differential equation is formed. Details of the proof can be seenin 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 11ATg /÷ATg + AT‘hdp÷dpATTE+EATdt dt2Figure 1.8: Normal Directrix Approach by DunwoodyThree constraints are necessary for a surface to be developable. They are stated asfollows using the above nomenclature:1. The normal directrix and generator vectors must be perpendicular.dpgDifferentiation with respect to t yieldsg2. The vector is of unit length.77=’.Differentiation with respect to t yieldsChapter 1. Introduction 12-;g3. The normal is invariant along a generator.( _;X g) - = 0Combining the constraints yields the following differential equation describing thenext consecutive generator:(;-‘\-g_16dt — 1dp dpdt dtIntegrating this differential equation using such integration routines as Runge-Kuttafourth 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 tobe of any practical use. Further constraining is necessary in order to control the behavoirof the developable surface.This approach is where the present research project started and expanded implementing new terminology and concepts which will be presented in detail.1.4 New Approach Research ObjectivesThe new approach research objectives were based on expanding the work initiated byDunwoody utilizing the normal directrix approach. This approach is unique in that itChapter 1. Introduction 13will always yield a solution. In this state, however, it is not very useful from a practicalengineering view point since it does not match two space curves. This leads to theresearch objectives presented in this thesis to further refine this technique.• The first research objective was to develop an algorithm in order to find a normaldirectrix 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 developablesurfaces and to generate the flat plate layouts and angles.• The final objective was to implement these algorithms using a modern computerlanguage and a popular CAD package in order to assess the practicality of theapproach.Chapter 2SplinesThis chapter is included because the material covered on splines contains the necessarybackground 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 particularcharacteristics which proved to be useful for a designer.2.1 Types of Splines ChosenWhen considering using a mathematical representation for space curves one can classifythem as of either parametric or non-parametric form. Non-parametric forms are usedextensively in various fields of mathematics and engineering. Non- parametric curves canbe further categorized as either explicit or implicit. The explicit non parametric form isusually expressed in the following form:y = f(x)where,x = independent variablef(x) = function of independent variabley = dependent variable14Chapter 2. Splirzes 15In this form multiple-valued or closed curves cannot be expressed. To overcome thisform, one usually uses the implicit form of the non parametric curve in the followingtypical form:f(x,y) —0where,f(x,y) = A function of both x and yOne typically calculates a point on the curve by calculating the roots of the equation.The approach can sometimes prove to be fairly computationally expensive. Implicitformulation is a very common form of non-parametric polynomials. Many formulationsused in engineering require higher than third order thereby making computations evenmore expensive when solving for roots of the equation.The non-parametric implicit formulation presents difficulties when being applied indefining such three dimensional objects as ship hull curves and surfaces. One typicalproblem that arises is when a vertical slope is encountered along the curve or surfaceresulting in an infinite numerical value. An infinite number cannot be used in a numericalcalculation when using a computer. Another problem arising in non parametric implicitform is that the positions are not distributed evenly along the curve or surface. Thisposes 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 representing closed curves and curves with multiple values of an independent variable. Parametriccurves are also axis independent. Parametric curves replace the use of geometric slopesChapter 2. Splines 16(which may be infinite) with parametric tangent and curvature vectors, which are neverinfinite. 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 polynomials 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 anddesired characteristics of the polynomials were investigated. Considering the types ofapplications and desired flexibility of the types and characteristics of curves desired, afairly exhaustive investigation of various types of splines was conducted. Special cubicpolynomials 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 behaviourof the splines yield cubic polynomials. Cubic polynomials are most often used becauselower-degree polynomials give too little flexibility in controlling the shape of the curve,and higher-degree polynomials can introduce unwanted “wiggles” and also require morecomputation. No lower-degree representation allows a curve segement to interpolate (passthrough) two specified end points with specified derivatives at each endpoint. Given acubic polynomial with its four coefficients,eg.f.’_ 43 j2x,1—a, +VX +c+four knowns are used to solve for the unknown coefficients. For example, the fourknowns might be the two endpoints and the derivatives at the endpoints. Other knownsmight be slopes or additional points[12]. It should also be noted that parametric cubicsChapter 2. Splines 17are the lowest-degree curves that are nonplanar in 3-D (three dimensions). You can seethis by recognizing that a second-order polynomial’s three coefficients can be completelyspecified by three points and that three points define a plane in which the polynomiallies. 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 usedin 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 developedBasis functions, and appropriate nomenclature on the various levels and types of curvecontinuities. No detailed analyses of the derivation of the splines used in this thesis willbe discussed in substantial detail since this work has already been done by such authorsas those previously cited. Only enough explanation of the nomenclature used in thisthesis to familiarize the reader with the concepts and characteristics of the various formsof the splines used will be discussed.Another point to mention is that local control of the 3-D curves was desired so thata curve segment is completely controlled by only four control vertices; therefore, a pointon 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:a1 a2 a3 a4 F_1b1 b2 b3 b4 F,P(i + t)= 3 i (2.1)C1 C2 C3 C4d1 d2 d3 d4Chapter 2. Splines 18where,P(t) = [x(t), y(t), z(t)j (2.2)andPoint = P(t) = [t3 t2 t i] [] [ (2.3)Tangent T(t) = = [3t2 2t 1 0] [ } [] (2.4)Curvature = C(t) = &t) = [6t 2 0 0] [1 [1 (2.5)2.2 Uniform Parametric, Geometric and Non-Rational Continuity RequirementsOne of the important properties discussed in such fields as finite elements and computeraided geometric design is of the mathematical techniques of shape representation. It istermed continuity. Continuity can be described as the highest level of differentiationwhich is continuous [3]. The types of continuity can be further categorized. Four typesof continuity are considered in this analysis. Each type of continuity will be explainedbriefly. The types of continuity chosen can be expanded but only two of what was thoughtto 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 usedin 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, thesetypes of splines appeared to be a likely choice. Non- uniform splines are not able to beexpressed in the format that was adopted in this reseach at this time.Another type of continuity requirement which was desired was parametric continuity [18]. Parametric continuity can be explain quite briefly with the aid of Figure 2.1. Ifthe nt1 derivative vector of two cubic curve segments are equal (ie. their direction andChapter 2. Splines 19magnitudes are equal) at the segments’ join point, the curve has nthdegree continuity inthe parameter t, or parametric continuity[12]. One would then state that if the directionand magnitude of [P(t)] through the th derivative are equal at the join point, thecurve is called C’ continuous. Figure 2.1 shows a curve segment S joined to three different curves with three different degrees of continuity ascertained by the superscript abovethe C. One should note that a parametric curve segment is itself everywhere continuous;the continuity of concern here is at the join points[12j.IC2x(t)Figure 2.1: First Three Levels of Parametric ContinuityIf two curve segments join together, the curve has GO geometric continuity. If thedirections (but not necessarily the magnitudes) of the two segments’ tangent vectors areequal at the joint point, point the curve has G’ geometric continuity. In computer- aideddesign of objects, G1 continuity between curve segments is often required. G’ continuitymeans that the geometric slopes of the segments are equal at the join point. For twotangent vectors TV1 and TV2 to have the same direction, it is necessary that one bea scalar multiple of the other. We then state the relationship that TV1 = k . TV2 withk> 0 [3j[12jy(t)Chapter 2. Splines 20One should note that in general, C’ continuity implies G’, but the converse is generally not true. G’ continuity is generally less restrictive than is C1, so curves can beG1 but not necessarily C1. However, visually, join points with C1 continuity will appearjust as smooth as those with C’ continuity as can be seen in Figure 2.2.[12j.In Figure 2.2 curve segments Q,, Q, Q join at the point P2 and are identical exceptfor their tangent vectors at P2. Q, and Q2 have equal tangent vectors , and hence areboth G’ and C’ continuous at P2. Q, and Q have tangent vectors in the same directionbut Q has twice the magnitude , so they are only G’ continuous at P2[12]Another type of continuity which one may desire is Rational Continuity [13]. Rationalcontinuity can be simply defined for “general rational cubic curve segments as ratios ofpolynomials:— X(t)— Y(t) — Z(t)2 6— 147(t)’— 147(t)’ z(— 147(t) . )where X(t), Y(t), Z(t) and W(t) are all cubic polynomial curves whose control pointsare defined in homogeneous coordinates. We can also think of the curve as existing inhomogenous 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 afourth 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. Forfurther reading one may refer to either Barsky and Hohmebar [13] or Foley [12jChapter 2. Splines 21y(t)ç,PFigure 2.2: Comparing Geometric and Parametric Continuity2.3 Uniform Non-Rational Tension Catmull-Rom SplineThe uniform non-rational Tension Catmull-Rom spline was chosen because it exhibitssevera’ useful features the designer may require [8] [91. First, it is an interpolatingspline, meaning that the curve passes through the points (control vertices). Second, itis in parametric form, meaning that it does not encounter singularities, only variationsof vector magnitudes. Third, it exhibits desired parametric and geometric continuityrequirements. Fourth, it has a global tension parameter which can further control theshape of the desired curve.The uniform non-rational Tension Catmull-Rom spline is easiest described in vectormatrix form. In this form it is a relatively simple task to imbed into a program. Thevector-matrix format is in a form in which not only the position along the curve canChapter 2. Splines 22be 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 splinesare presented.This pseudo-standard form shows that the spline has a local influence of the controlvertices as any chosen position. At any one particular location along the curve the controlvertices 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 particularposition along the curve the only infuence is from P_1,P, P:+iandPj+2.Taking into account the end conditions of the splines also had to be considered. Theapproach taken was to create Phantom points. Given the control vertices vector below,Phantom end conditions were formulated. Detailed analysis of the derivation can bereferred to in the appendices D.1 D.2. The vectors located below show the indexing ofthe control vertices and how end conditions are treated (phantom points).Pi-1PiPi+1Pi+2For the initial condition P0 the control vertices vector has the following form D.l:Chapter 2. Splines 23PiPiPi+1i+2For the end condition P,_1 the control vertices vector is in the following form D.2:PiPiFi+1Pi+1The following page shows the vector-matrix form of the uniform non-rational TensionCatmull-Rom spline. The spline is parameterized with respect to the parametric variablet and has a global tension variable, 6.—2.0/3 4.0— 2.0/3 2.0/3— 4.0 2.0/31 4.0/3 2.0/3 — 6.0 6.0 — 4.0/3 —2.0/3P(t) = t3 t2 t’ 1—2.0/3 0.0 2.0/3 0.00.0 2.0 0.0 0.0Chapter 2. Splines 24Pi-1PiPi+1Fi+2= Tension parameterTo give the reader a good comparison as to the behaviour of the Tension CatmullRom spline the following figures are included. The first figure, figure 2.3(a) shows thespline with a tension value of 0.5 shown relative to interconnected line segments. Forthe tension Catmull-Rom spline with a value of 0.5 this defaults to a traditional cardinalspline.The next figure, Figure 2.3(b) shows little difference but is relaxing the spline tensionwhen given a tension value of 1.0.Finally the Tension Catmull-Rom spline in Figure 2.3(c) shows how it nearly contoursto the inter-connected line segments shown in the figure. If the tension parameter is givena value of 0.0 then it becomes line segments.2.4 Uniform Non-Rational Beta-SplineThe 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 controlvertices. 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 doesChapter 2. Splines 25(a) Tension at 0.5 Catmull-Rom (b) Tension at 1.0 Catmull-Rom(c) Tension at 0.1 Catmull-RomFigure 2.3: Catmull-Rom SplinesChapter 2. Splines 26not have any singularities occur. Third, it also retains desired parametric and geometriccontinuity requirements. And, fourth, it has global bias and tension parameters whichfurther 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) t2 t’ 1]—2.0/3? 2.0(132 + i3? + /3? + j3) —2.0(132+/3? +/31 + 1.0) 2.01 6.0/3? —3.0(132 + 2.0/3? + 2.0/3?) 3.0(132 + 2.0/3?) 0.0—6.0/3? 6.0(/3?— /3k) 6.0/3k 0.02.0/3? (132 + 4.0,3? + 4.0/3i) 2.0 0.0Pi-1PiPi+1= i3 + 2.0/3? + 4.0/3? + 4.O/3 + 2.0= Bias/32 = TensionChapter 2. Splines 27Like the Tension Catmull-Rom spline the Beta-spline exhibits the same localizedcontrol vertices influence. Again, if one were to choose a particular position along thecurve the closest control vertex would only be influenced by the preceding one and thenext two, eg. [P_1,P, P+i, P÷21.Taking into account the end conditions of the spline also had to be considered inthe same fashion as the Catmull-Rom spline. The same approach was taken to createthe Phantom points. Given the control vertices vector below, Phantom end conditionswere formulated. Detailed analysis of the derivation can be referred to in the appendices 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+i6—2.Of3FiPi+1The end condition control vertices vector for the Beta-spline at P(n-l) is:Pi-1PiPi+12.013?P +(4.0+4.0/31+/32)Pt-i-i3—2.0Chapter 2. Splines 28The first of the Beta-spline figures shown below reveals the spline compared to interconnected line segments. The first figure, figure 2.4(a) shows the curve near the controlvertices, exhibiting the characteristics of an approximating spline. The bias is at 1.0 andthe tension is at 0.0. With these values the Beta-spline degenerates to a B-spline, whichexhibits 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 tensionat 25.0. At these values one can see how this spline can also resemble inter-connectedline 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 anegative 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 isleft at a value such as 0.0 the curve exhibits the following behaviourThese figures shown above try to give the reader some idea of the capabilities of thesetype of splines and how one can use the features each spline possesses. These are but afew of the types of splines which are now being developed. Each of these types of splineshas features which the reader should be aware of in order to maximize the benefits theyhave 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 BetaFigure 2.4: Beta SplinesChapter 3Creation of a Normal Directrix from Two Space CurvesThe normal directrix approach will always yield a smooth developable surface in thevicinity of the normal directrix, so long as the normal directrix is itself smooth. Unfortunately, it is not always clear what shape the normal directrix should take in order thatthe resulting surface should meet the requirements of the designer. With respect to therequirements of the designer, the definition of a developable surface from two edges issuperior. The objective of the present work is to create an algorithm which will shape adevelopable surface to lie close to a pair of edge curves. The developable surface will bespecified in terms of a normal directrix in order to ensure that the surface is smooth. Itis not expected that the developable surface will contain the two edge curves, only thatit 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 onepoint on the first generator, then constructing a curve which lies within the surfaceand is perpendicular to all generators. In addition, extra construction tools, in theform of differential equations, were created in order to better control the normaldirectrix solution.2. Once a normal directrix has been computed, it can presumably be smoothed out toyield a smoother developable surface without departing greatly from the originalequation.30Chapter 3. Creation of a Normal Directrix from Two Space Curves 313. Nolan’s approach of matching the cross products between the generator and thetangents to each of the edge curves can be expressed in terms of a differentialequation.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, thenthe values for those n points can be computed to yield a spline which lies close tothe normal directrix derived from the differential equations.The approach taken to try to match a developable surface to two edge curves resulted in formulae modelled by differential equations. The differential equation versionof the conventional method, named the modified conventional approach, created by Dun-woody and Konesky was used as an initial guess in order to utilize additional differentialequations to solve for directrix control vertices.3.1 Modern Approach Utilizing a Single Normal DirectrixA normal directrix can be computed for any developable surface by starting at one pointon the first generator, then constructing a curve which lies within the surface and isperpendicular to all generators. This very powerful technique was created by Dunwoodywhich is termed in this thesis as the modern approach. In addition, extra constructiontools, in the form of differential equations, were created in order to better control thenormal directrix solution.This modern approach, given a normal directrix and a start generator position, willalways force a developable surface to be created. This solution is underconstrained,however, and further refinement was deemed necessary in order to better control theChapter 3. Creation of a Normal Directrix from Two Space Curves 32behaviour of the function.The formulation is in vector differential equation form and uses the ODE class integratioll routine which implements Runge-Kutta 4th order. Constraints for this developablesurface differential equation are given in the next subsection. For a more thorough analysis one can refer to Appendix B for the derivation.3.1.1 Constraints Governing Modern ApproachThe differential equation developed by Dunwoody [10] termed the modern approach involves three constraints which define a developable surface. They are as follows:1. The generator must be of unit lengthie.2. The vector normal is invariant along a generatorie.(gx).=03. Vectors must be perpendicularie.dg dp_ d2pdt dt dt2 gThe three constraints can be seen in Figure 3.1 below in vector form.Using the three constraints which define a developable surface yields the differentialequation in simplified form:d (. “dg— dt2 g, Pdt — \dt dt)Chapter 3. Creation of a Normal Directrix from Two Space Curves 33g‘1Illi:dt dt dt2Figure 3.1: Derivation of Developable SurfaceThis equation works fine as is but further constraining is necessary in order to makethis solution practical. For example, given the theory presented so far we can generate adevelopable 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 desiredfinal generator position. We now move to the next stage in the development of alignmentof the end generators.3.1.2 Alignment of End GeneratorsOnce a normal directrix has been computed, it can presumably be smoothed out to yielda smoother developable surface without departing greatly from the original equation. Inaddition, a much more realistic and desireable condition is where the user gives a normaldirectrix, a start generator position and a final generator position and the configurationadjusts itself to conform to the newly added constraints.Chapter 3. Creation of a Normal Directrix from Two Space Curves 34Change in Angle With Respect to Unit Motion of a Control VertexThe problem was approached by looking at the problem using the perspective of observingthe location of the directrix, the position of the generator and separate the vectors into in-plane and out-of-plane components. Another differential equation was formulated whichaccumulated information of the rate-of-rotation of a generator with respect to motionof 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 correctaccordingly to match the end generator. Figure 3.2 shows the relation of the originaldirectrix and the corrected directrix as well as the amount of change that is necessary.— ___. ———OriginalDirectrixModifiedDirectrixFigure 3.2: Vector Locations and Corresponding AnglesA detailed analysis of the derivation of the differential equation can be referenced inAppendix C if the reader wishes to look further. A brief summary of the highlights of thederivation is needed here in order to familiarize the reader with new concepts which arebeing introduced with the theory and the nomenclature of the user controllable designvariables.The angle Theta, herein referred to as 9, is the out-of-plane rotation of G’ withrespect to G. The angle Phi, is herin referred to as 4’, and is the in-plane rotation of G’N’TChapter 3. Creation of a Normal Directrix from Two Space Curves 35with respect to G. The variable “a” is a parameter describing a directrix control vertexcomponent.The rate of rotation of a generator with respect to changes in one of the parametersof the directrix curve is written as:dadtThe desired expression is the rotation of a generator with respect to changes in oneof the parameters of the directrix, which is written as:dOdaThe desired expression can be defined in terms we can derive, namely:= -•N (3.2)where, N is the unit normal defined as: (3.3)N = Txg (3.4)but, (3.5)ddO— d(dgN 36dtda — dt”darewriting gives, (3.7)d28 — d2g dN38dadt— dadt dt dtUsing numerous identities and proofs, the user can look at the derivation in detail inAppendix C, the resulting differential equation which contains both the in and out-ofplalle components result in the following formulation:—39dadt— (“) . )Chapter 3. Creation of a Normal Directrix from Two Space Curves 36Using several identities and constraints that have already been presented the finalform simplifies down to the following relation:O(x) &pdadt — (2 (3.10)\dt dt’The Alignment Process and the Concept of MobilityThe 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}2 d2p1N—14!dt= 1N—1 ( X ) dadt{5} dt (3 11Jo dadt Jo\dt dt’dadt{5+1.O}dadt {6+2.O}where is defined as:dadta bcdefgh= 3.0t2 2.Ot 1.0 0.0 (3.12)i j k 1 6+1.0m n o w 6+2.0Chapter 3. Creation of a Normal Directrix from Two Space Curves 37The alignment process may take several passes, ie. from 0 to N-i, a correction, thenfrom 0 to N-i, etc.The procedure will be described and the successive passes usually reduce the error byone decimal place per pass (ie. iteration).After the first accumulation of information along the spline from t = 0 to t = N1, a few constraints are added. The first desired constraint of the movement of thecontrol vertices is that the end control vertices can only move in the direction along eachcorresponding end generator. The second from the end control vertices are constrainedto move both in the direction of the end generators and in the direction of the startand end tangent vectors respectively. Both ends are calculated in the same way, so for abetter 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 thedirection of the starting generator. The end condition constraints are also influenced bythe 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.IPantomPointFigure 3.3: Control Vertices and End Phantom PointAt the end of the spline the end two control vertices of the summed terms have toChapter 3. Creation of a Normal Directrix from Two Space Curves 38be modified because of the phantom end condition constraints. The end control verticesare influenced by the following condition which describes only the B-spline criterion:2F — P1. From this relation we have to add again to itself for the P2 control vertexof information and subtract of the P2+’ component from itself. This is in essence thetechnique which is applied to the end conditions of the control vertices for each typeof spline being used.The concept of mobility is quite simple and provides a further constraining on thebehaviour of the alignment process. Each control vertex of are assigned a mobilityvalue, which is a constant. A mobility value of 1.0 means that the correspondingvertex has full mobility and that it is free to adjust that control vertex as governed bythe equation. At the other extreme a mobility value of 0.0 indicates that the vertexis not to be moved during each iteration of the alignment process.3.1.3 Analysis of ResultsMany tests were done with the modern approach and then it was incorporated as afoundation to build upon. The approach was found to be very robust but specificallycontrolling its behaviour became the dominant problem. It was also noted to have onemajor 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 solutionexisted already, thereby, no need of the modern approach would be required.Mobius Strip Demonstrating FlexibilityAs an example which demonstrates both the flexibility of the modern approach ie. thecontrol vertices were not permanently fixed in position, (NB. only the end conditionswere 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 39if one was to follow an edge it would continue infinitely along the surface travelling alongboth sides. For clarity, the reader should refer to Figure 3.4 in this chapter.Figure 3.4: Robustness of Modern Approach Showing Mobius StripIn order to create the developable mobius strip the very end control vertices werefixed. As mentioned earlier a new concept was introducted. This concept was termedmobility and is explained in ‘the next section below with figures included to help clarifythe explanation.Controlling Alignment with MobilityThe concept of mobility was developed as a first tool to control the behaviour of themodern approach. It can simply be defined as local weighting of the change in anglewith respect to unit motion of each control vertex. Each control vertex has informationrecorded about it by the function described in equation 3.10. When a mobility value of1.0, unity, is assigned the class entity which retains the information about a particularvertex 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 canuse a range of values for each vertex ranging from 0.0 to 1.0 if we wish to “fine tune” orChapter 3. Creation of a Normal Directrix from Two Space Curves 40more accurately control the behaviour of the iterative solving procedure in order to alignwith the end generator.The example of the developable mobius strip shown previously and now shown intwo views in figure 3.5 was constructed with 6 control vertices. Each control vertexwas assigned a mobility. The first and last vertex were given a mobility of 0.0, ie. nomovement allowed, and the other four were given full mobility of 1.0, ie. full freedom tobe corrected.(b) Developable mobius strip view 2Figure 3.5: Mobius StripThe developable mobius strip was one of the first examples where mobility was necessary to create a specific type of shape. The same 6 control vertices used to construct themobius strip were all initially given mobility of 1.0, and run to see what type of solutionwould result. The resulting figure 3.6, demonstrates what type of solution results when(a) Developable mobius strip view 1Chapter 3. Creation of a Normal Directrix from Two Space Curves 41no constraining is implemented.Figure 3.6: Modern Approach Showing Full Mobility of All PointsOne should note that figure 3.6 is a valid solution. The end generators do align withthe same line that passes through both the start and end generator positions. Clearlyone 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 onesolution could result; an entire family of solutions is possible with the single directrixapproach unless further constraining is included.Problems arising when Surface has Small In and Out-of-plane CurvatureIn the modern approach a problem arises when the surface has small in and out-of-planecurvature , ie. the surface is almost flat. This problem is located where the alignmentprocess takes place. When the accumulated values of the change in angle with respectto motion of the control vertices is very near zero, a subsequent calculation involvesChapter 3. Creation of a Normal Directrix from Two Space Curves 42dividing 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. Ifthe plate is very close to being flat the computer being used to do the analysis will yielda floating point error due to division by zero.A tolerance or threshold value was implemented which would make the resultingfloating point operation equal to zero if the total in and out-of-plane curvature was lessthan 1.0 x 10_b. This was deemed necessary in order to improve the robustness of thealgorithm. This results in the modern approach to accomodate when the surface is verynear flat when initialized. A sample is shown in figure 3.7 where the surface is flat anda corresponding correct solution results.HFigure 3.7: Provision Made for When Surface Very Near Flat3.2 Constraints Defining Modified Conventional Approach (Modified Nolan’sApproach) in Order to Create a Normal DirectrixNolan’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.Nolan’s approach can give good approximations of developable surfaces. This approachexpressed in differential equation form and with additional constraints could be used asan approximation for creating a close fitting normal directrix.Before a normal directrix can be computed, it must have information as to whereChapter 3. Creation of a Normal Directrix from Two Space Curves 43it should be located relative to two space curves that contain a developable surface.Information in the form of modelling the Conventional Approach by differential equationslead to the formulation of Normal Directrix Control Vertices. This material is presentedwith the Modern Approach Utilizing a Single Normal Directrix because it is used to derivethe normal directrix control vertices.In figure 3.8 the two outer space curves, referred to as f(u) and g(v), are used indefining the modified conventional approach. In this approach the first two of the threeconstraints stated by Nolan [17] are exactly the same. The third definition is also usedbut is modified to include an additional term is formulated in order to determine wherealong the generator the normal directrix should lie. These steps are now explained inmore detail.NiCuef(u)Cueg(v)Figure 3.8: Relating Modified Conventional Approach Space Curves and Normal DirectrixTangentPlaneIn figure 3.8 the first two constraints are graphically shown stating that each spaceChapter 3. Creation of a Normal Directrix from Two Space Curves 44curve must have each tangent vector perpendicular to its normal. Also, the cross-productsbetween 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:dg— X gendvFigure 3.9: In-plane and out-of-plane componentsThe third relationship is derived from the out-of-plane components and the in-planecomponents leading to their next respective locations along the space curves. The thirdrelationship is stated in equation 3.16 and is explained here below.dt—(genxni)duwhere tand ii;;= —Xgendu= (g(v)—f(u))=(g—f)(3.13)(3.14)(3.15)Out ofd2—— dt— dg(gen X n2) .(3.16)Chapter 3. Creation of a Normal Directrix from Two Space Curves 45The numerator dot product relation in equation 3.16 represents the out-of-plane component. The denominator relation represents the in-plane component. Only magnitudesare needed since direction is constrained to being along each space curve’s tangent andnormal vectors. At a particular position along the space curve the out-of-plane directionis located along the normal vector and curvature vector. For example, on the space curvef(u), i and are in the directions of out-of-plane curvature. For in-plane curvatureone 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 magnitudeswould be simplest to equate with respect to the next parametric calculation along thespace curve since the directions are constraints in the formulations. These relations arethen equated to two space curves. By equating these relations to two space curves aknown developable set of generators (or rulings) can be first approximated. The equating 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 thatequation 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 thespace 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 conventional 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) (3.17)From the two space curves, f(u) and g(v), the generator is found as follows:Chapter 3. Creation of a Normal Directrix from Two Space Curves 46= (g—f) (3.18)Equation 3.17 is then differentiated with respect to t in order to establish a differentialequation for the normal directrix intersection with the generator. Differentiating yields:dp dfdu dgdv da= (3.19)Using the definition that, . = 0 and using this on equatioll 3.19 results in thefollowing:dp df du dg dv__da= (3.20)The normals of each space curve is also calculated using the following relationships:= df (3.21)= x gen (3.22)Relating the out-of-plane curvature with the in-plane curvature of the two spacecurves as mentioned previously we get:_;: = . (3.23)(gen x ni). (gen x n2).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 47dv— 2n (1 nt du 3 24dt—2ri,,dtRearranging these relationships to solve for , and ij gives the following relation:—1d2fdu—+(genxni).—325dt— 2n 2n (. )dv2—-(genxn2).—dvdv— 2n ndu326dt— t ndt—.--- du dg —s dv= (l—a)—.gnj+a--.gen..(3.27)gen gen3.2.1 Offset of Normal DirectrixOne feature which was found useful for the designer to have as a variable to adjust isthe initial positional offset of the normal directrix. This relation is shown again here asequation 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 theuser may be able to further fine tune the surface to desired specifications. If this is notof any concern it defaults to the parametric value, 0.5, which is the midpoint betweenthe two space curves.Chapter 3. Creation of a Normal Directrix from Two Space Curves 48f(u)g(v)Figure 3.10: Positional Offset of Normal Directrix3.2.2 Allowment of Different Number of Control VerticesDuring the initial design stages one criterion was noted, namely that the number ofcontrol vertices for each space curve and the normal directrix may be desired to bedifferent from each other. This was noted and a further relationship was establishedwhich allows for different control vertices for any of the curves. The relationship ispresented 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 491.0 du 1.0 dv 1.0=— (3.29)Equation 3.29 states that each parametric differental variable, and , contributesone-half times its total number of control vertices. Ignoring the total number of controlvertices for every curve results in the equation being written as:+ 1.0 (3.30)The resulting relationship for the space curves is as follows:= (331)dt flj ndtFurther reading detailing the derivation can be referred to in G.3.2.3 Present Problems with Current Solution and Improvisation ImplementedTests including equation 3.16 proved promising. One problem noted, is shown in Figure 3.11 below when the surface is close to becoming flat. A phenomenon which wetermed “fanning” occurs where the out-of- plane component becomes very small relativeto the in-plane-component yielding undesirable characteristics. This is shown below inFigure 3.11. One can see from Figure 3.11 that intuitively if a plate is close to becomingflat, the generators should lie nearly perpendicular to the tangent vectors of each spacecurve, resulting in square flat plates for a flat surface. This problem was noted and addressed by incorporating an out-of-plane tolerance which would default to parallel lineswhen the surface fell below the user-specified tolerance. This is a temporary solution andfurther research is being directed to address this problem in future work.Chapter 3. Creation of a Normal Directrix from Two Space Curves 50Figure 3.11: Fanning Occurring When Developable Surface Nearly FlatFor a more detailed example showing where this phenomenon occurs in practice seethe examples in Chapter 7, entitled “Demonstration Examples”. The three examplescited were a simple conical developable, an arctic series fishing vessel, and the UBCseries fishing vessel.3.3 Utilizing Modified Conventional Approach to Approximate a Single DirectrixIf the normal directrix is to be described by a spline with n control points, then thevalues for those n points can be computed to yield a spline which lies close to the normaldirectrix 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 numberof control vertices for the normal directrix is still to be determined and may vary. Thenumber of desired normal directrix control vertices is specified by the variableFrom the above mentioned relationships another differential equation is needed inorder to solve for the location of the normal directrix control vertices. The relationshipwas found to be in the form of finding the minimum for an equation relating the normaldirectrix in terms of an error. This relation relates p(t), the known spline, and r(t), thespline control vertices we wish to solve for, is as follows:N—ierror = Ip(t)—r(t)I2dt (3.32)Chapter 3. Creation of a Normal Directrix from Two Space Curves 51We then differentiate equation 3.32 with respect to the control vertices and solve forthe relation to determine a minimum. This is shown as follows:Oerror= 2JN1 3- (p(t) — r(t)) dl (3.33)N—2=YJ.‘(p(j+s)-r(j+s))dS (3.34)0= 0 (3.35)Details of the analysis can be seen in Appendix F. The initial relation showing theequation is as follows:T8ij—1 Cj_1=[A]s2 [ . 1 ] [A] dS (3.36)6ij+1 S cj+18ij+2 1 C3+2Cj_1 6ij—1C2 N—2 —1 N—2 1 ‘5ii= [ (ii)] r(j+s) [3 2 s 1 ] [A] d3.37)Cj+l 5ij+1C,2 5ij+2Chapter 3. Creation of a Normal Directrix from Two Space Curves 523.3.1 Equations Yielding Approximated Normal DirectrixCollecting equations 3.25, 3.26, 3.27, equation 3.37 and presenting them below showshow all of the variables are sequentially coupled and integrated using the class ODE(Ordinary Differential Equation) solver.______ —1d2du—+(genxn1).—3 38dt— 2n 2n . )dv2(genxn2).—dvdv—2-i,, ndu339dt— nt ndt-du g dvcia — (1 — a)— gen + a— gen (3 40)d1Cj_1 sij—1C, N—2 —1 N—2 1= [(ii)] yj r(j+s) [3 2 [A] d3.41)Cj+1 sij+145ij+2In the computer program all of the terms are collected as the ODE solver integratesfrom 0.0 to N-L The algorithm tested concluded to be quite robust. No singularitieswould occur. There are still, a few problems that have to be addressed which still affectthe solution.Chapter 3. Creation of a Normal Directrix from Two Space Curves 533.3.2 Results and Present ProblemsCombining the equations above and integrating using the ODE class to generate the generators for the developable surface and the control vertices for the normal directrix yieldsthe developable conical-like shape shown as an example in Figure 3.12. The directrixshown in the figure has actually two directrices. Both happen to coincide exactly for thisexample. One directrix follows exactly as the surface is being integrated. The otherdirectrix results from solving for the normal directrix control vertices.Figure 3.12: Conic-Like Section with Flatness Tolerance Specified at 0.001The next figure shows the same conical-like shape without a flatness tolerance specified. The problem of fanning is apparent at both ends. Note how the two directrices aredifferent. The one which looks invariant is the control vertices which have been solvedfor 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 relatesthe apparent tangent with the current generators.The following figures display the modified conventional approach used to approximateChapter 3. Creation of a Normal Directrix from Two Space Curves 54Figure 3.13: Conic-Like Section With No Flatness Tolerance Specifiedthe normal directrix control vertices which is then used by the modern approach. Theexamples 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 currentconfigurations.In Figure 3.15 both versions appear to give better results. The modified approach hasthe tolerance set higher , 0.01, and displays the fanning problem near the stern beingreduced.In figure 3.16 is an overlayed closer comparison of the two methods with the currentconfigurations.In Figure 3.17 both versions appear to give even better results. The modified approachhas the tolerance set higher , 0.1, and displays the fanning problem near the stern isalmost eliminated.In figure 3.18 an overlayed closer comparison of the two methods with the currentconfigurations is shown.Chapter 3. Creation of a Normal Directrix from Two Space Curves 55(a) UBC Series Surface (main deck chine andchine 1) body plan(b) UBC Series Surface (main deck chine andchine 2) profileFigure 3.14: Developable Surface Procedure Comparison, tolerance 0.001Chapter 3. Creation of a Normal Directrix from Two Space Curves 56(a) UBC Series Surface (main deck chine andchine 1) modified approach body plan(c) UBC Series Surface (main deck chine andchine 1) non-optimized body pian(b) UBC Series Surface (main deck chine andchine 2) modified approach profile(d) UBC Series Surface (main deck chine andchine 2) non-optimized profileFigure 3.15: Developable Surface success, tolerance 0.01Chapter 3. Creation of a Normal Directrix from Two Space Curves 57(a) UBC Series Surface (main deck chine andchine 1) body pian(b) UBC Series Surface (main deck chine andchine 2) profileFigure 3.16: Developable Surface Procedure Comparison, tolerance 0.01Chapter 3. Creation of a Normal Directrix from Two Space Curves 58(a) UBC Series Surface (main deck chine andchine 1) modified approach body plan(b) UBC Series Surface (main deck chine andchine 2) modified approach profile(d) UBC Series Surface (main deck chine andchine 2) non-optimized profile(c) UBC Series Surface (main deck chine andchine 1) non-optimized body planFigure 3.17: Developable Surface success, tolerance 0.1Chapter 3. Creation of a Normal Directrix from Two Space Curves 59(a) UBC Series Surface (main deck chine andchine 1) body pian(b) UBC Series Surface (main deck chine andchine 2) profileFigure 3.18: Developable Surface Procedure Comparison, tolerance 0.1Chapter 4Optimization of a Normal Directrix4.1 ObjectiveOnce the modern approach was deemed fairly stable the next criterion were then approached and listed as follows:• To better match a developable surface defined by a normal directrix to a pair ofspace curves.• To create an optimization function which is the mean square distance from eachspace 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 integratedmean square distance.4.2 Optimization FunctionCalculation of the normal distance between a point on a space curve and the developable surface, see Figure 4.1, can be obtained by three definitive relationships. Thiswill 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 product60Chapter 4. Optimization of a Normal DirectrixFigure 4.1: Relating the Distance Between Space Curves and Surfacedpof the tangent,--and the generator, g(t).-dpn(t)=xg(t)61• 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 theclosest point from the surface to a point in space would be a specified distance, 1,macecurvesDeJves tosurfacedevelopablesurfacedlrectrixChapter 4. Optimization of a Normal Directrix 62normal 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 bedifferentiated with respect to the independent variable t in order to create a differentialequation which can solved.Differentiating yields the following:drdu —ds —dl dp dg dn—g(t)—n(t) = (4.1)We also have two other known relations which can be substituted, namely:-fU__.\— 1dp (42)dtdi dt= (4.3)From these three relations we can solve for three equations and three unknowns, withthe unknowns being the following dependent differentials:• The differential dependent parametric position along a space curvedudtThe second form involves using the three definitive relations and continuing to findrelations 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- nTaking the dot product with vector g:( Jdu ds_dl_g) — g • g — gThis simplifies to the following:dsdi= dt di di= g+sg+lg( du(4.5)(4.6)(4.7)(4.8)(4.9)Chapter 4. Optimization of a Normal Directrix 63• The differential scalar dependent position along the generator g(t)dsdi• The differential scalar dependent position along the normal n(t)dldiWe then solve for the following relation:4E-i--cit dt. dt dt-‘-i (4.4)cit dty dty dtydt dtz dtz dtzdrdux —g —ndrduy —gy flydrduz —g —nChapter 4. Optimization of a Normal Directrix 64Similarly, dot product with n, (4.10)and simplifying yields: (4.11)dl (dr ..idu—= I—n I— (4.12)dt du dtSimilarly, dot product withy (4.13)(E - 414du dt) dt — dt dt S dt dt di diSimplifying yields: (4.15)i . — -. (s g + 1 n)du_________________________di ( 4.16du di4.2.1 Integral Chosen to Minimize and Simplex Parameters UsedThe key parameter chosen to calculate, the distance from a point on a space curve tothe closest point on a surface, has to be further related to a function which describes theentire surface. Calculation of a least squares approximation for the area under the curveswas chosen. One can refer to figure 4.2 for a visual representation of the area under thecurves. 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 formwould 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)b= I Iv(t)Idt (4.18)Ja&= f v(t)dtJad ptJav(r)dr=v(t)d pt= dtfa Iv(r)Idr = Iv(t)IdS = v(t)dtArea= J F2 Idh1 I 2 d1] dt[1j+l2Chapter 4. Optimization of a Normal Directrix 65AlecUndO.xvedevelopablesurfacedlrecfrlxFigure 4.2: Area calculated under space curvesdSdtLet h1 and h2 be arc lengths forThe resulting least squares type(4.19)(4.20)(4.21)(4.22)the two space curves r1 and r2 respectively.approximation for the total area would be:(4.23)Chapter 4. Optimization of a Normal Directrix 66From equation 4.23 the area between the space curves and the developable surface 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 offreedom for the non-linear optimization method model. The non-linear optimizationmethod used is the Downhill Simplex method. This is discussed later in this chapter.By solving all the differential equations’ independent and dependent variables fromthe start of the developable surface to its end, the non-linear optimization function canattempt to minimize equation 4.23.Another key point is that the dependent and independent variables sum to give tendegrees of freedom. The nonlinear optimization method also needs ten different solutionsto start with. This problem was solved by moving one control vertex degree of freedomfrom the directrix by unit length. From these initialization procedures the DownhillSimplex method was able to start solving these equations.4.2.2 Present ProblemsBoth approaches stated above describing solving for the dependent variables will produce solutions for surfaces which do not have major changes in the rate of change ofthe 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, andchanging in direction, the modern approach may have difficulty in the alignment processwhen trying to conform to the surface. Cases may arise, as shown in figure 4.3, where acorresponding normal on the surface which should match a position on the space curvewill not occur. This may be due to the constraint that both the space curves and thesurface are of finite length and the corresponding normal on the surface will point toan out-of-bounds position. One other possibility is that a localized change in curvatureorients itself midway during the iteration process and another normal of infinite lengthChapter 4. Optimization of a Normal Directrix 67occurs.The first method of matching the modern approach to two space curves can also failat a prior stage when calculating the inverse of the matrix. Reasons for a non - solutionare generally for the same reasons.DIsonce LI frornspoceletonormalFigure 4.3: Possible inherent failure due to certain geometry4.3 Downhill Simplex MethodMany non-linear optimization methods as cited by Jacoby [14] would be possible toimplement, but, a more suitable and generally more stable and robust method, theSimplex method [16], was chosen. A major reason for this preference was the descriptioncited in Press [20], which briefly described its multidimensional capablities.Chapter 4. Optimization of a Normal Directrix4.3.1 Explanation of The Downhill Simplex Method68A brief description of the Downhill Simplex Method in Multidimensions is presented hereto 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 MethodThe downhill simplex method has three parameters which define the moving behaviourof the simplex. These parameters are mentioned here because changes in their valuesChapter 4. Optimization of a Normal Directrix 69may be necessary. The three parameters are defined by Nelder [16] as a, 6, and.These parameters correspond to: reflection, contraction and expansion. The reflectionparameter, a, is a positive constant, which is a scalar multiplication constant whichmirrors the point through the simplex. The contraction parameter, 3, is a constantwhich values lie between 0 and 1. It is a ratio of the distance of the point relative tothe simplex centroid. The expansion coefficient,,is greater than unity and it is a ratioof the current point to the centroid with a point along the line joining the point to thecentroid. For a more detailed explanation refer to Nelder [16]. For detailed coded form ofthe Downhill Simplex Method, programming and tests performed, refer to Appendix L.4.3.2 Results and Present ProblemsAs stated previously, surfaces which resemble close to developable surfaces to begin withwill be more likely to have a solution that closely matches the original space curves. Beloware two examples which demonstrate both simple and difficult solutions encountered.4.4 ExamplesA few examples are shown in this section to help explain the various phenomenon of theresults of the theory created this far in the research.4.4.1 Example of Single Surface with Conical PropertiesThe first example demonstrates how the non-linear optimization can successfully solve afree-form surface with conical properties. Shown below in figure 4.5 are two views of aconical like surface.Note that due to the symmetry of the object in figure 4.5 the generators also followsymmetrically along the surface.Chapter 4. Optimization of a Normal Directrix 70(a) Conical-like developable surface view 1 (b) Conical-like developable surface view 2Figure 4.5: Conical-like developable surface4.4.2 Example of UBC Series Demonstrating Present ProblemsThe UBC Series vessel was posed as the immediate goal to make into a developable shipseries. In figure 4.6 the surface tested was using the 1st and 2nd chine as the spacecurves.As seen in figure 4.6 the optimized version had difficulty with the adjustments due tothe 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. Themodified approach has the tolerance set low , 0.001, and displays the fanning problemnear the stern.Chapter 4. Optimization of a Normal Directrix 71(a) UBC Series Surface (Chine 1 and 2) non-optimized(b) UBC Series Surface (Chine 1 and 2) optimized(c) UBC Series Surface (Chine 1 and 2) comparisonFigure 4.6: Developable Surface FailureChapter 4. Optimization of a Normal Directrix 72(a) UBC Series Surface (main and chine 1)modified approach body plan(main and chine 1)(b) UBC Series Surface (main and chine 2)modified approach profile(d) UBC Series Surface (main and chine 2)optimized profile(c) UBC Series Surfaceoptimized body planFigure 4.7: Developable Surface success, tolerance 0.001Chapter 5Intersection of Developable Surfaces and Flat Plate LayoutOnce developable surfaces have been created, it is necessary to define the edges of thosesurfaces by intersection with adjacent surfaces.An additional consideration is once the bounded developable surface has been createda flat plate layout is necessary. A derivation of the solution is also presented here.5.1 Intersection of Developable SurfacesOnce the developable surfaces have been created there is the requirement that they mustintersect without gaps or spaces between them. Shown below is the vector representationfor initial conditions of the surfaces and their orientation relative to each other. Infigure 5.1 we see the various vectors that are dependent upon each other. We can showthe vectors dependence on each other as follows:Figure 5.1 presents the independent parametric variable I. The two dependent parametric variables are u and v. The independent space curve is p(t); one dependent spacecurve which is a function of u is f(u); the other dependent space curve is g(v), which isa function of v. The curves p(t), f(u) and g(v) are normal directrices representing thedevelopable surface of interest, and the adjacent surfaces on each side. The intersectioncurves joining two developable surfaces are r1 which is relating f(u) and p(t) and theother intersection curve r2 relates g(v) and p(t).A detailed derivation of the equation used to calculate the intersection of the developable surfaces are located in Appendix I73Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 74Agft:472p(\X*/ ...•S b(k dhh(v)Figure 5.1: Intersection of Three Developable Surfaces5.1.1 Derived Equations Yielding Intersections of SurfacesDerivation of the equations used to calculate the intersection of the developable surfacesare as follows:r1= Pt + it (5.1)dr1 dp dg ds1—k- = +sit--+--gt (5.2)Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 75•. +(Edtgt—-dt —___dg gt’dpdp d)-- SS4. Igu Idfudu — df dfdrdu-dr1duQfudgu du-;u- = + S2 — + fJBut,= f + s2ugu, ands2ugu =— fand, r1 = p+sgSubstituting gives,dfdfdu(dhdhd2fu du=dt _S2u)(5.3)(5.4)(5.5)(5.6)(5.7)(5.8)(5.9)(5.10)(5.11)(5.12)(5.13)(5.14)Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 76/ —s __..(dpdu_+it+-du—— (5.15)dt ‘—1df df d2f—.——-..(pt+sitgt_fu))-- ---i(5.16)dt — du dtr2 = Pt + S2tgt (5.17)- -dr2 — dp dgt ds2--- —+ +—— gt (5.18)—k /_ .... /____.dp (d1i dgt (dh!__•ã._xv) +S2rcit — /___. (5.19)(dIig Ija x g)I____. \dg YtIdp (5.20)•Id2pdtE\ (5.21)r2 = h + (5.22)Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 77dr2 — dr2v523dt— dvdtdr2= (5.24)r2 h +s1 (5.25)r2—(5.26)r2 p + sä (5.27)dr2 dh dv (dh dh d2h= (5.28)(Idt 2tdt dvdv J(dh d1z d2h—1i)dgdgdv 530dtdvdt5.1.2 Results and Present ProblemsThe intersection of developable surfaces was tested and proved to be fairly stable whenwell behaved surfaces were created.Conical Type Surface TestedAn example shown in figure 5.2 displays very tightly intersecting surface without anyproblems. 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 78Criterion of Phantom SurfacesThe criterion of phantom surfaces was introduced to provide a secure mechanism forpreventing interconnecting surfaces from having “cracks” or “leaks”. The conical-typesurfaces presented in figure 5.2 show surfaces which have been created using the modifiedapproach. These all appear to be developable but only the inner two have been confirmed as developable using the modern approach. The two outside surfaces were used asphantom 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 approachand yield aligned surfaces with the intersection algorithms. The cracks will be moreevident 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 surfacesICI(a) Conical-type surfaces view point (1,1,1)\\\\\\\(b) Conical-type surfaces front viewFigure 5.2: Conical-type Surfaces IntersectionsChapter 5. Intersection of Developable Surfaces and Flat Plate Layout 805.2 Flat Plate LayoutIn order to convey the information of a developable surface into a useful form, the fiat-plate-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 AnglesOnce a developable surface has been created, the shape and position of the consecutivesegments that make up the developable surface need to be known. Below is the derivationof the algorithm used to provide the information of the shape of the segments in the xy plane and the angle and rate of change of bending of the segments relative to eachprevious one. For a more detailed explanation of the derivation one should refer toAppendix J/!2 LhtdtdtFigure 5.4: Flat plate layout derivationThe resulting differential equation is used to find the new tangent and correspondingChapter 5. Intersection of Developable Surfaces and Flat Plate Layout 81(a) 3-D view of developable mobius strip (b) Corresponding flat plate layout(5.31)Figure 5.5: Developable surface and flat plate layoutThe flat plate layout corresponding to the numbered plates in figure 5.5 has datalisted in table 5.1.point on the x-y plane where the next point of the plate is located:12a q \dt2 dt) dt UP=5.2.2 Format of OutputThe example of the flat plate layout for the developable mobius strip can be seen infigure 5.5Chapter 5. Intersection of Developable Surfaces and Flat Plate Layout 82Plate no. Parametric value out of plane angle (deg)0 0.00 4.641389 0.9751381 0.25 13.924170 8.4518002 0.50 23.206951 19.5084273 0.75 32.489716 27.8896694 1.00 32.115269 21.8432465 1.25 29.928911 13.1609476 1.50 38.130260 19.5311477 1.75 56.255440 60.3593798 2.00 55.513470 41.1225559 2.25 34.989170 8.16899710 2.50 21.341158 3.40450111 2.75 14.529543 2.68255812 3.00 12.068573 6.00717213 3.25 9.672798 13.69153914 3.50 10.510071 9.68131415 3.75 15.940498 9.44085816 4.00 18.051476 19.69595317 4.25 12.895594 21.22638718 4.50 7.738149 11.44155219 4.75 2.579468 4.39382520 5.00 -0.016314 0.398113Table 5.1: Flat plate dataChapter 6Choice of Computer Language and CAD ProgramWhen this thesis work was initiated, one of the objectives was to initially create generictools in a very portable modular language and then apply them to a well establishedopen ended CAD package that is used throughout industry. At the present stage of thework both objectives seem to have been met.6.1 Computer Language ChosenThe choice of which computer language to program was based on a number of designconstraints. The language chosen should be flexible in case the method of design had tobe altered. The three designs considered were:• Top-down structured design• Data-driven design• Object-oriented designOf these three types of design, data-driven design was initially dropped because of thenature of the design tool. Because of the financial constraints on the project and thesoftware languages available on the platform chosen, the top-down structured design wasinitiated. This was then taylored into a modular top-down structured design. Afterthe release of a popularly supported compiler for the platform was chosen, the designwas then switched to an object-oriented design. The language initially chosen was ANSIC, best supported at the initial time of development by Microsoft for the PC platform.83Chapter 6. Choice of Computer Language and CAD Program 84Later, after an initial attempt with Walter Bright’s, Zortech C++, the switch was madeto Borland’s C++, which proved to have excellent support tools, a good developmentenvironment and better technical support.6.1.1 OOP - Object Oriented ProgrammingWhen it was realized that a true portable object-oriented language was available, thedesign was taylored using influences from OOP, object oriented programming, OOD,object oriented design, and OOA, object oriented analysis. First a definition of each isin order:• OOP Object-Oriented ProgrammingObject-oriented programming is a method of implementation in which programsare organized as cooperative collections of objects, each of which represents aninstance of some class, and whose classes are all members of a hierarchy of classesunited via inheritance relationships [5].• OOD Object-Oriented DesignObject-oriented design is a method of design encompassing the process of objectoriented decomposition and a notation for depicting both logical and physical aswell as static and dynamic models of the system under design [5].• OOA Object-Oriented AnalysisObject-oriented analysis is a method of analysis that examines requirements fromthe perspective of the classes and objects found in the vocabulary of the problemdomain [5)Chapter 6. Choice of Computer Language and CAD Program 856.1.2 Portability to Different PlatformsOne major consideration when doing software development for an application is thestandardization and portability of the code being developed. Up until a few years ago, oneof the most standardized and portable languages was ANSI C. Now with more featuresand better type checking, ANSI C++ appears to be the successor.Porting some of the sample computer code mentioned below to different compilersand 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 theSUN SPARC UNIX platform. No major problems were encountered.6.1.3 C++: Classes and FeaturesThe language chosen to best implement the work being done was the C++ language. Thislanguage proved to be very much more than just a superset to the C language. Many ofthe class features, when used in its entirety, better shape the programming environmentinto a true object oriented program and design. A little explanation is presented in thenext sections as to how the tools were used to better inform the reader as to the approachtaken.Class vector*iaclude <iostream.h>*ifndef _VECTOR_*define _VECTOR_en direction { I, Y, Z };class matrix;Chapter 6. Choice of Computer Language and CAD Program 86class vector {mt n;static jut default_length;float *eleent;public:vector0 {ndefault _length; ele.entnev float En ; };vector(int din) {ndi;eleentnew float[n];};vector(int di., float x);vector( coast vectorl vvector0;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 operatorfriend vector operator*(const matrixk a, coast vectork b);vector& operator*(float a);vectort operator*(const vector& a); II e1e.entby—element nltiplicationfriend vector operator/(const vectork v, float a);friend istrea.& operator>>(istreamk, vector&);friend ostrea.k operator<<(ostreamk, const vector&);*endifThe class vector is the name of a class which enables the user to use N dimensionalvectors for calculations. A C++ feature, termed operator overloading enables the userto write equations almost exactly as one would write them by hand.Chapter 6. Choice of Computer Language and CAD Program 87Included below is the vector.h class header file which shows how the class is initiallydesigned. Some explanation is in order at this stage to clarify to the reader these basicfeatures which are used throughout the development of other classes used in developingthe 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. Below is an optional enum direction { X, Y, Z }; which enables the user to declare X as a0, Yas land Zas 2.The next line class matrix; declares the matrix class before class vector.h can bedeclared, 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 classname. 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 areprivate: by default.The next line, whose member access control is public: shows some very interestingfeatures of C++. The first constructor, depicted by vector() (n=3; element=new double[n]);, which shows a default size, n, for the class to allocate memory. vector() is anoperation which is empty of a type. This case is where the type is hidden and somemechanism must be provided for a user to initialize variables of that type. Located onthe next line is, vector(int dim) n=dim; element=new double[n], which is a constructor oftype double. Moving down one line, vector(const vector & v), is a copy constructor. “Acopy constructor for a class vector is a constructor that can be called to copy an objectof class vector that is, one that can be called with a single argument of type vector. Acopy constructor is generated only if no copy constructor is declared”. The next line,iector() delete element;;, is contains what termed a destructor. A destructor automatically deallocates dynamic memory when a class goes out of scope. Moving to the nextChapter 6. Choice of Computer Language and CAD Program 88line contains one possible case of operator overloading. The line, doubleF_4 operator[](inti)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 inthe following form, float operator[](mt i) const{returri element[i];};. This is necessarysince the compiler complains when several operations are done successively and invoketemporary storage variables. For example:vector x,y,z,a,b,c,d;.x,y,z,a,b,c,d . . .get assigned valuesChapter 6. Choice of Computer Language and CAD Program 89thenx (y*z) *a* (b*c) *d;cout<<x;The above list of code shows how the operator jj is used as well as the << streamoperator is overloaded to standard output for the class vector. Syntactically this code ismuch simpler and far more legible to the reader. This code is simple and easy to followthereby 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 bya member function within classvector. So, for example, in order to find the size ii of adeclared class of type vector, one has to invoke either friend mt size (vector éY v)returnv.n;; or 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 couldencounter 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 ispassed is left for the reader to investigate. The second operator, vectors operator=(constmatrix& a);, equates two objects of class type vector and matrix. This enables an equating of classes of different types. The next example operator,friend vector operator *(constvectors a, const vectors b);, requires a friend to be declared. This code is optional sincethe same feature could be invoked by the following, vectors operator*(vector a);.Chapter 6. Choice of Computer Language and CAD Program 90Class matrixIn 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 N1JLL);friend matrix identity(int);friend matrix inverse(coast matrix& a, float* detNULL);friend istreamk operator>>(istreama, matrixi);Chapter 6. Choice of Computer Language and CAD Program 91friend ostreant operator<<(ostreaak, const satrixa);friend hit operator<(const .atrix&, conat •atrixt);friend zatrix abs(const aatrix&);I;*endifOne can see the similarity between class vector and class matrix and how they caninteract.Class Curve and Surfacetinclude “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 = new vector[n];};Curve(Curvek a);CurveO{delete[n] point;};void realoc(int na);vectork operator El (hit i){return point Ii) ;};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);Chapter 6. Choice of Computer Language and CAD Program 92static void initbasis(char sptype’2’, double b11.0, double b20.0);static atrix BEB;tendiftifndef_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:Chapter 6. Choice of Computer Language and CAD Program 93ofstrea. file;Filelisto *next;ofstrea.k operator[](iut n);class Filelisti{public:ifstream file;Filelisti *next;ifstream& operatorlj (jut a);*endifThe following classes, Curve arid Surface also implement vector and Curve. Thesesort of implementations were able to make the code more legible and apply the testingof 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 Initiallyvoid rk(double tO, double& del_t, double ti, vectorl x,vector (sf)(double t, vectort x), vectork err_scale);#endifThe above code class ode, ordinary differential equation solver, also shows anotheruseful combination of class building to contribute to a library of useful class tools.6.2 Computer Platform SelectedChapter 6. Choice of Computer Language and CAD Program 956.2.1 Pc 386/486 EnvironmentThe development of the theory for this thesis and the validation via progranmiing wasimplemented on a PC. The PC, personal computer, platform was decided upon becauseit 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 hardwarecan be interfaced. With PCs becoming faster and less expensive every year most of thespeed constraints of the thesis programs were relaxed. At the present state of the thesisresearçth and design, the programs execute fast enough to become interactive. Commercialviability should be considered.6.2.2 computer Language Selection in this PlatformAt the start of this thesis the computer language chosen was the C programming language.This language proved to be sufficient for initial implementations. Further developmentof computer languages, such as C++, proved to offer much more sophisticated featuresand modularity of code that would offer more efficient development cycles. This provedcorrect and the theory, code, test, development cycle time was greatly reduced.Many other features of the C+++ language also proved invaluable towards implementing the theory. One of the most used features of the C++ language was the useof operator overloading. This enabled vector calculus equations to be written in code inalmost the same convention one would use when writing by hand. Other features whichwere discussed above all contributed to more legible code and as a collection of usefultools which other applications would ensue.Chapter 6. Choice of Computer Language and CAD Program 966.3 CAD Program SelectedDuring the course of development of the theory a decision was made as to whether ornot writing a computer platform independent CAD package should be undertaken. Afterassessing the existing CAD packages used on several computer platforms, the conclusionwas that one would be wiser to develop the special code to an existing CAD package. TheCAD package chosen was assessed as the most popular and efficient to develop upon. ThisCAD package was AutoCAD. During the course of development of the theory, ongoingassessment was done on AutoCAD’s market position and development tools. This provedto be benefitial to this project as well as it reflected upon us how strong a commercialproduct this was in the CAD market as well as their desire to move to other hardwareplatforms.AutoCAD was reviewed and concluded to be the most popular CAD package in industry and contained the best development design tools out of all of the CAD packagessurveyed. Initial development started with AutoCAD release 10 on the 386 personal computer. 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 andsolve most of our earlier problems. Now pending publishing of this thesis ongoing development is being done towards commercialization of developable surfaces as a third partydeveloper package. Much effort is needed in order to create user interactive tools sophisticated enough to ensure useful interactive CAD tools that the user would find intuitiveto use.6.3.1 CAD Program Environment and Open ArchitectureThe CAD program environment using AutoCAD ADS development tools proved to bewhat was needed for developable surfaces. The ADS development tools are C languageChapter 6. Choice of Computer Language and CAD Program 97calls which are registered and loaded as applications for AutoCAD. The ADS development tools have proved to be of high quality and as an open architecture to build andexpand upon the existing tools contained in AutoCAD. This open architecture has enabled AutoCAD features and newer releases to grow very rapidly. The writing of ADSdevelopment tools in ANSI C language has also enabled AutoCAD to be ported to otherplatforms very quickly and easily with minimal effort when compared to other languages.6.3.2 Present LimitationsSo far the present limitations of AutoCAD and ADS development tools are few. Some awkwardness has occured when implemented AutoCAD ADS applications tools with some ofthe C++ modules. True C++ OOD, OOP techniques can be affected by some of theC 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++ classesand methods. Existing design changes are being implemented quite easily into the AutoCAD ADS environment. Autodesk has been very supportive and provided encouragementfor our ongoing research and development implementing C++ to ADS.Chapter 7Demonstration ExamplesIn this chapter there are four example sections which display the benefits and problemsstill facing the present design. The first section describes the power and flexibility of thebase module which, when unconstrained, can create more difficult geometric anomaliesin computational geometry.The other sections show more realistic applications and the difficulties that arisetrying to “best fit” two or more space curves to the exact geometry desired. One mustbe made aware that from any two or more space curves, a developable surface may notexist; only a closest fitting one may be a solution. The reader should be made aware thatthis is a design tool, and it should also provide the user with feedback as to what typeof geometries are potentially developable.7.1 Developable Mobius StripThe 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 problemwhen working from a directrix and two generators only.This geometric anomaly is interesting to view and cite as an example and perhaps asmathematical art, but provides no potential commercial benefit to the industrial designer.The next few sections exemplify more practical examples in commercial areas such asmanufacturing. When “harder” constraints are involved new problems arise. Examplescould have been cited where these design algorithms would prove to the reader that they98Chapter 7. Demonstration Examples 99(a) Developable mobius strip view 1 (b) Developable mobius strip view 2Figure 7.1: Developable mobius StripChapter 7. Demonstration Examples 100work very well, but, they would also prove to be misleading. This thesis was presentedin a very objective engineering manner which shows both the positive and negative sidesinherent in design.7.2 Simple Conical DevelopableA simple conical developable surface was created in order to verify and check compliancewith known developable surfaces, fanning and intersection of surfaces.Shown above 7.2 are surfaces which intersect and verify results.7.3 Arctic Fishing Vessel7.4 UBC Series Fishing VesselReferring to figure 7.5 the ubc series was found to be the most difficult to model since thebow contains compound curvature hull form in certain regions from station 0 to station2.0.Chapter 7. Demonstration Examples 101(c) Conical-type surfaces view point (1,1,1) (d) Conical-type surfaces front view(a) Conical-type surfaces view point (1,1,1) (b) Conical-type surfaces front viewFigure 7.2: Conical-type Surfaces IntersectionsChapter 7. Demonstration Examples 102(a) ARCTIC vessel view point (1,0,0)(c) ARCTIC vessel view point (0,0,1)(b) ARCTIC vessel view point (0,4,0)(d) ARCTIC vessel view point (-1,-1,0.2)Figure 7.3: Arctic Vessel Conventional ApproachChapter 7. Demonstration Examples 103(a) ARCTIC vessel view point (1,0,0) (b) ARCTIC vessel view point (0,-1,0)J// ‘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 ApproachChapter 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 ApproachChapter 8Conclusions and Recommendations8.1 ConclusionsThe research objectives were• to develop an algorithm in order to find a normal directrix such that the resultingdevelopable surface lay close to two space curves representing desired edges of thedevelopable surface.• to create an algorithm to intersect developable surfaces and to generate the flatplate layouts and angles.• to implement these algorithms using a modern computer language and a popularCAD package in order to assess the practicality of the approach.Normal Directrix from Two Space CurvesAn algorithm was created to compute a normal directrix for a developable surface froma pair of space curves. Differential equations were derived for defining a set of generatorslying 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 approximatedby a parametric spline.To ensure that the rulings between the ends of the space curves were also generatorsof the developable surface, the tangents at the ends of the space curves were alignedso that the cross-products between the end rulings and the two tangents at each end105Chapter 8. Conclusions and Recommendations 106were co-linear. The control verticies of the normal directrix were adjusted so that thedevelopable surface started from one end of the normal directrix aligned with the oppositeend generator.The algorithm created yielded undesireable fluctuations in generator direction whenthe surface was nearly flat. A threshold limit for the curvature was used, below whichthe generator direction was held constant. Practical examples with hard-chine ship hullsdemonstrated the robustness of the approach.A non-linear optimization technique, the downhill simplex method, was implementedto further refine the shape of the developable surface. Limited success was achieved.Intersection of Developable Surfaces and Flat Plate LayoutAn algorithm was derived to define the edges of a developable surface by intersectionwith adjacent developable surfaces. The method proved to be robust, with the exceptionwhen two adjacent developable surfaces are nearly co-planar.Differential equations were also derived for creating a flat plate layout of a developablesurface, including axes of bending and plate curvatures.ImplementationThe C++ programming language was used to implement the algorithms. This enabled anobject-oriented and modular approach, for example 3-D space curves were implementedas a class of objects with member functions such as position, tangent and curvature asa function of the independent parameter. As a class of objects, details of the implementation of curves were transparent to or isolated from the rest of the program. Otherfeatures of this programming language enabled the main program segments to be writtenclearly and concisely in a mathematical format.Chapter 8. Conclusions and Recommendations 107AutoCAD was selected as the host CAD package because of its wide acceptance,multiple platforms and development tools. Only preliminary steps have been taken todate 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 adjacentsurfaces when nearly co-planar.• Implement non-uniform rational beta spline and non-uniform rational tension CatmullRom 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 ComputerGraphics. In Proceedings at the Pacific Northwest Section of the Society of NavalArchitects and Marine Engineers, pages 1 — 21, January 1987.[3] B.A. Barsky. 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 BoundaryValue Problems. John Wiley & Sons, 1977.[7] J.C. Clemens. A Computer System to Derive Developable Hull Surfaces and Tablesof Offsets. Marine Technology, 18(3):227— 233, July 1981.[8] T.D. DeRose and B.A. Barsky. Geometric Continuity and Shape Parametersfor 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 Geometric Constructions for Catmull-Rom Splines. ACM Transactions on Graphics,7(1):1— 41, January 1988.108Bibliography 109[10] A. B. Dunwoody. Computer Aided Design of Developable Surfaces. not yet published, 10, 1989.[11] editor E.V. Lewis. Principles of Naval Architecture, 2nd ed. Volume 1,2,3, Societyof Naval Architects and Marine Engineers, 1989.[12] J.D. Foley, A. VanDam, Steven K. Feiner, and John F. Hughes. Computer GraphicsPrinciples 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. PrenticeHall, 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. ComputerJournal, 7:308— 313, 1965.[17] T.J. Nolan. Computer-Aided Design of Developable Hull Surfaces. Marine Technology, 233 — 242, April 1971.[18] J.C. Beatty R.H. Bartels and B.A. Barsky. An Introduction to Splines for Use inComputer Graphics and Geometric Design. Morgan Kaufmann Publishers, 1987.[19] D. F. Rogers. B-Spline Curves and Surfaces for Ship Hull Definition. Society ofNaval Architects and Marine Engineers, 1977.[20] S.A. Teukolsky W.H. Press, B.P. Flannery and W.T. Vetterling. Numerical Recipesin C. Cambridge University Press, 1988.Appendix AMathematical Notation for Partial DifferentiationNotation 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 hasthree terms:dp — a7 dx+op dy+dzdt Oxdt 9ydt OzdtThe use of the straight “ d” denotes derivatives of only one variable.For example:The Derivative means F’(t)• op awhile means f1(x,y,z) = —f(x,y,z)rIOAppendix BDerivation of Developable SurfaceB.1 Constraints which Define the Developable SurfaceThe definitions presented below were from previous work and on going research..i!LtTg g+.JTgdt dt dt2Figure B.1: Derivation of Developable SurfaceVectors Must Be Perpendiculard d d(g+ T)(+LT) 0= 0)dg dp(+g = 0)dg dp—cPpdt dt — dt2 g111Appendix B. Derivation of Developable Surface 112Vector g is of Constant LengthdgThe Normal is Invariant Along a Generatorbax(r,t) g(t)rp(t)tFigure B.2: Derivation showing normal is invariant along a generatorx(r, t) = p(t) + rg(t) where r is a scalardx dp dg=dxb=-;--=gardp dgn = normal=x=(+r)xgdp dg= (-jjxg)+(r--xg)Appendix B. Derivation of Developable Surface 113dp dgBut (- x g) x g) ze. are paralleltherefore ( x g). = 0B.2 Proof Using Constraintsdg_ (.g)dpdt (2)dtLet B be the constant we need to satisfy:dg— Bdt dtCONSTRAINTS1. Vectors must be of unit lentgth.g•=O==g•B=02. Normal is invariant along a generator.(gx). 0=(gx).B = 0Identity: A. (B x C) = B. (C x A)B5.(gx) = 0Bg.(x) = 0Appendix B. Derivation of Developable Surface 1143. Vectors must be Perpendicular.dg dp— d2pdt dt —B — d2p“dt dt’— dt2 g(2.— _dt2 g-\dt dttherefore: = — g) dpdt (.)dtAppendix CDerivation of Rate of Rotation of Generator Differential EquationThis derivation also includes some new terminology.A differential equation is derived which calculates the rate of rotation of a generatorwith respect to changes in one of the parameters of the directrix curve.NN’ThetaPhiDirectrixFigure C.1: Vector locations and corresponding anglesThe angle Theta, herein referred to as 0, is the out of plane rotation of G’ with respecttoG.The angle Phi, is herein referred to as , and is the in plane rotation of G’ withrespect 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 parameters115Appendix C. Derivation of Rate of Rotation of Generator Differential Equation 116of the directrix curve is written as:d20dadtThe desired expression is the rotatior of a generator with respect to changes in oneof the parameters of the directrix, which is written as:dOdaThis desired expression can be defined as:da dawhere, N is the unit normal defined as:N Txgtherefore,£4--(.Ndtda — dt\dagiving,d28 — d2g dg dNdadt — dadt dt dlFrom the derivation of a developable surface:\dt2dlVdt dtwhere T is of unit length:T = dt/.4RVdt dttherefore,dadt da da /Vdt dt Vdt dtAppendix C. Derivation of Rate of Rotation of Generator Differential Equation 117The second term in the above equation disappears because the derivative of a constantis zero.therefore,— (g)dTdadtVdt dt(C.2)dadt \da Jdt dtOne of the next terms to define is the first term of the second part of= çT + ONdawhere,dT=_—— gaacombining gives,— ( . g) T + ON (C.3)One of the next terms to define is the second term of the second part ofUsing the differentiation product rule on one of the definitions:N=Txgyields,dN /dT ‘\ / dg’\—- =x g) + x (C.4)Appendix C. Derivation of Rate of Rotation of Generator Differential Equation 118Combining Equation C.3 and Equation C4 to form the second part of equation C.1,(—(.g)T+oN). ((xg)+(Tx ))= (_ . g) (T. (4 x g)) + ((_z . g) (T. (T x +9(N.(xg))+6(N.(Tx))several terms drop out of Equation C.5:0 (N. (Txff))= BT, hence, T x TdtAnother term also drops out:o(N.(x))Using Identity A. (B x C)(4JZ(g xN)),butTresulting in -- Tsince, T remainsOne other term drops out:(—.g) (T.(Tx))butdtThereforedg dNda dtUsing a vector triple products dot= 0 since,= 0= B.(CxA)= gxNunit0,length= BT, resulting in (T x T) = 0and/dT ‘\( IdT——g) T.—jj-xgcross identityda dt —(C.5)Appendix C. Derivation of Rate of Rotation of Generator Differential Equation 119IdT ‘\ IdT= -_.g)-.(gXT)using the fact that N = T x gfdT ‘\ IdT= da dtThen,d20— ( g) (dT N + (dT (dT Ndadt—da ) da gj dt-___da — da’ 1daVdt dtUsing the quotient rule for derivatives:dT V di di \. dadt) dt \. dt dt) dadt dida\dt di/ ,12d2______dadt dt dp=dt dt dt dt)- t’_(dadtdt).Nda /4.4z (dtThe second term drops out because N = 0ThereforeUIN \dadtda—Vdt didl’ d1N•— = N..1 didtVdt dtbut Using quotient ruleAppendix C. Derivation of Rate of Rotation of Generator Differential Equation(2\dt di\cit di(-— dadt dt2)\dt dt— ‘dt2 ) dadt\di di120—V dt di dt2) di \. dt di) dt2 dtIV (dd\dt dtd2— dt2Vdt dt(.& /— \di2 di)(E.& dt\di di)Again, the second term drops(ç.N)Vdt didpout because — N = 0dt(.dadi g/2.42Vdt dtdTNdtdTNdtdTdadTad20dadtd20dadt(.& /dadi dt) dp—______(42 at‘di di)-‘-.N1•g)— i..dt2 J ‘dadi (\dadi g di2(.\dadiVdi diCombining terms yields,(42.2’ (1\di di) di diUsing Identity: (A. C)(B. D) — (A. D)(B. C)= (A x B) . (C x D) yields,( X • (N x g)(x)gxNBut, A•(BxC)=(CxA).B— \.di2 di— 3(&di di)d2pdadtAppendix C. Derivation of Rate of Rotation of Generator Differential Equation 121d29______— dt2 dt)—___. (C.6)dt dt)Appendix DTension Catmull-Rom SplineThe Catmull-Rom Spline matrix with a tension parameter, /3, is shown below. Thephantom point end conditions are shown on the following pages.—2.0/3 4.0 — 2.0/3 2.0/3— 4.0 2.0/31 4.0/3 2.0/3 — 6.0 6.0 — 4.0/3 —2.0/3P(t) = u3 u2 u1 1—2.03 0.0 2.0/3 0.00.0 2.0 0.0 0.0Pi-’PiPi+1/3 = Tension parameter122Appendix D. Tension Catmull-Rom Spline 123D.1 Phantom point, P(O), at 0—2.0/3 4.0 — 2.0/3 2.0/3— 4.0 2.0/31 4.0/3 2.0/3 — 6.0 6.0 — 4.0/3 —2.0/3P(t) = 0.0 0.0 0.0 1.0—2.0/3 0.0 2.0/3 0.00.0 2.0 0.0 0.0Pi-1Pi= Fl_iFi+iFl+2P1-P10.0 1.0 0.0 0.0 = Fl__iPi+iFl+2Appendix D. Tension Catmull-Rom Spline 124—2.0/3 4.0 — 2.0/3 2.03 — 4.0 2.0/31 4.0/3 2.0/3 — 6.0 6.0 — 4.0/3 —2.0/3F(0) = 0.0 0.0 0.0 1.0—2.0/3 0.0 2.0/3 0.00.0 2.0 0.0 0.0FiPiFi+1i+2Appendix D. Tension Catmull-Rom Spline 125D.2 Phantom point, P(1), at n—2.0/3 4.0— 2.0/3 2.0/3 — 4.0 2.0/31 4.0/3 2.0/3 — 6.0 6.0 — 4.0/3 —2.0/3P(t) = 1.0 1.0 1.0 1.0—2.0/3 0.0 2.0/3 0.00.0 2.0 0.0 0.0Pi-’PiPi+1Pi+2Pi-’Pi0.0 0.0 1.0 0.0 = Fi+2Pi+1Appendix D. Tension Catmull-Rom Spline 126—2.0/3 4.0— 2.0/9 2.03— 4.0 2.091 4.0/3 2.0/3 — 6.0 6.0 — 4.0/9 —2.0/3P(1) = 1.0 1.0 1.0 1.0—2.0/3 0.0 2.0/3 0.00.0 2.0 0.0 0.0Pi-1P1Pi+1Pi+1Appendix EBeta-SplineP(t)= [ 2 1 1 ]—2.0/3? 2.0(132 +/? +i3? + /9i) —2.0(132 +i3? +/3i + 1.0) 2.06.0/3? —3.0(132 + 2.0/3? + 2.0/3?) 3.0(132 + 2.0/3?) 0.0(5—6.0/3? 6.0(/3?—6.0/3k 0.02.0/3? (/32 + 4.0/3? + 4.O/3) 2.0 0.0P1-P1Pi+1= 132 +2.0/3?+4.0/3?+4.0/3i+2.0= Bias/32 Tension127Appendix E. Beta-Spline 128E.1 Phantom point, P(O), at 0P(t) = [0.0 0.0 0.0 1.0]—2.O/3 2.0(/32 + /3 + i? + i3i) —2.0(,62 + /? + /3 + 1.0) 2.06.0/9—3.0(/92 + 2.09 + 2.0/3?) 3.0(132 + 2.0/3?) 0.0S—6.0/3? 6.0(13?—/3) 6.03 0.02.0,6? (/32 + 4.0/3? + 4.0/3k) 2.0 0.0Pi-1Pi=Pi-1Pi+1Fi+2S =/32+2.0/?+4.0,B?+4.0i91+2.040Appendix E. Beta-Spline 129Pi-11.0 Pig 2.0/3 (/32 + 4.0/3 + 4.0/3k) 2.0 0.0 =i+1Pi+2— (/32 + 4.O/9 + 4.O/31)P + 2.0P÷1(6—2.O,8)Appendix E. Beta-Spline 130P(0) = [0.0 0.0 0.0 1.0]—2.O/3 2.0(132 + i? + I? + /3k) —2.0C82 + /? + th + 1.0) 2.01 6.0/3? —3.0(132 + 2.0/3? + 2.0/3?) 3.0(82 + 2.0/3?) 0.0S—6.0,3? 6.0(13? — /3) 6.0,3 0.02.0/3? (/32 + 4.0/3? + 4.Ofli) 2.0 0.0(I32-I-4.OI3?+4.o31)P+2.oPiS—2.O3?PiPi+1Appendix E. Beta-Spline 131E.2 Phantom point, P(1), at nP(t) = [1.0 1.0 1.0 1.0]—2.0,6 2.0(/32 + /3 + /? + 3) —2.0(132 + /? + i3 + 1.0) 2.01 6.0/3 —3.0(132 + 2.0/3 + 2.0/3?) 3.0(/32 + 2.0/3?) 0.0S—6.0/3 6.0(/3— /3i) 6.0/3k 0.02.03 (/32 + 4.0/3? + 4.0/3k) 2.0 0.0Pi-1Pi= Pi+2Fi+1Fi+2S = /92+2.0/3+4.0/3?+4.0th+2.0Appendix E. Beta-Spline 132Pi-’1.0 Pi0.0 2.0/3 (4.0/3? + 4.0/3k + /32) 2.0 = Pi+2Fi+1i+2— 2.0/3F, + (4.0/3? + 4.06 + f32)P11i+2 8—2.0Appendix E. Beta-Spline 133P(1) = [1.0 1.0 1.0 1.0]—2.0/3? 2.0(/32 + i? + i? + /3) —2.0(/32 + /? + /3 + 1.0) 2.06.0/3? —3.0(/32 + 2.0/3? + 2.0/3?) 3.0(/32 + 2.0/3?) 0.045—6.0/3? 6.0(/3? — /3k) 6.0/3k 0.02.0/3? (/32 + 4.0/3? + 4.03) 2.0 0.0Pi-1P1Fi+12.Of3?1’ + (4.O+4.O31+P2 )P116—2.0Appendix FDerivation of Normal Directrix Control VerticesIt was found necessary to solve for a desired number of normal directrix control verticesfrom two space curves. The two space curves must be the same type of splines as thenormal directrix in order to minimize the complexity of calculations. The number ofcontrol vertices of each space curve can be different from each other as well as bothnumbers 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 involvingp(t), our known spline, r(t), the spline we wish to solve for, in the following error equation:error= 1N-1Ip(t) — r(t)12 dt (F.1)where we wish to minimize the integral and solve for the Contro Vertices, C:öerror= 2JN1- (p(t)— r(t)) dt (F.2)N—2 ld (+s)=dC, (p(j+s)—r(j+s))dS (F.3)= 0 (F.4)J’ dp(j+ s) pj + s)dS=2j’dp(j+s). r(j + s)dS (F.5)134Appendix F. Derivation of Normal Directrix Control Vertices 135where the the above terms are defined below in vector form:Ci-’Cip(j + s) = 3 2.[A] (F.6)Ci+1Ci+26ij—1ã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 = ATBT expanding for triple products (ABC)T =CT (AB)T = CTBTAT we get the following relation:Tãerror=[A]T [3 2C2dS (F.8)s1 C2Appendix F. Derivation of Normal Directrix Control Vertices 1361N—2= f r(j + s) [3 2 s 1 ] [A] dS (F.9)0‘5ij+1ij+2 jSimplifying some of the terms yields the following formTF j_’l 1 r 1Co Iôerror N2 1 [3 2 S 1 ] dS[A] (F.lO)ac = I [A]f3=01 II 5ij+l I I [ ICN1]Sij+2j [1T‘5ijl— 76541 r 1N_i 1 L 43I I IIc0 Iiii I I1 1 1 1 I5432 I= I I [A]T I [A] I j (F.11)j=0 . I61j+2 1 1 1 1432N—2= (i,j)1’ (F.12)j=0Appendix F. Derivation of Normal Directrix Control Vertices 137Cowhere 1’ (F.13)CN1N—2= (i,j) [1’] (F.14)j=o5ij—1N—2 1r(j+s) s 1] [Al dS (F.15)5ij+18ij+2(F.16)‘5i3—1N—2 —1 N—2 1[Fj = [(ii)] j r(j+s) [3 2 1] [Al dS (F.17)sj+16ij+2Once the control vertices have been solved for, the desired end condition constraintsmust 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 138where P are space curve control vertices and C are solved for control vertices.Also,(P1- F0) (F.20)= IF’—P01(FM — PN) (F.21)gN= FM—FoIAppendix GModified Conventional Approach DerivationA few trials of the normal directrix method yielded results which were very dependentupon the initial position of the normal directrix when trying to match the normal directrixmethod developable surface to two space curves. Further testing revealed that if a goodinitial guess of a normal directrix could be achieved to match to two space curves, thenlittle optimization or correction is necessary and a closest developable surface could befound.Below is a figure relating the variables used in the derivation of solving for an initialnormal directrix control vertices given the control vertices of two space curves.I.t)g(v)Figure G.1: Orientation of Space Curves and Directrix139Appendix G. Modified Conventional Approach Derivation 140From Figure G. 1 f(u) and g(v) refers to the space curves and p(t) refers to the normaldirectrix. 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)dfdu= (l_a)Tdgdv .da+ a--- + gen(G.2)(G.3)From Figure G. 1 we also calculate the normals at these points, namely:nl = ---Xgenandg= j—Xgen(G.4)(G.5)Now, relating the out-of-plane curvature with the in-plane curvature of the two spacecurves we get the following:-;dtX n1). —duAlso, relating the number of control vertices between the three curves and howand are related we have:Rearranging, gives:1.0dv — 2n 10dt—nt du— 2n dt (G.8)= I 211 Idv(gen x n2).—dvi(G.6)1.0 du2n dt1.0 dv+-(G.7)Appendix G. Modified Conventional Approach Derivation 141Rearranging, Equation G.6 to solve for and and gives:____. —1d2fi• IIdu t nt (genxni).1 IIS(G.9)d2g IIfl2-IIdg I(genxn2).— )dv’ /dv — 2n,, ndu (G.1O)-df d dgda (1—a).gen+aj---.gendt (G.11)Appendix HRelating the Space Curves to the Developable SurfaceTwo forms were developed:1. The first sets up the relations and solves for the equations2. The second calculates the parameters directlyThe First ApproachThere are primarily three constraints which determine the distance from a point onthe developable surface to the closest point on a space curve.They are as follows:n(t) dp(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 oft, which is the key parameterized variable. In order to solve for these we try to solve theproblem in terms of a differential equation relating the variables u,s,l, to t. This is done42Appendix H. Relating the Space Curves to the Developable Surface 143by partial differentiation of the function r(u). The differential equation is with respectto t is as follows:drdu dp dg ds dn dl=Rearranging into a useful form,drdu ds dl dp dg dn-g--n =There are only three unknowns to solve for given the initial conditions relating eachspace curve to the developable surface:1) At t = 0.0, 11 = 0.0 and 12 = 0.02) At t = 0.0, ul = 0.0 and u2 = 0.03) At t = 0.0, q(0.0,sl) rl(0.0) and q(0.O,s2) = r2(0.0)4) si = (q(0.0, si)—p1(0.0)) g(t)s2 = (q(0.0, s2)—p2(0.0)) g(t)and,and can be evaluated,dg - (Y)dpdt-\dt dtJdn — 1d2p 1dp dg— Xgj+XThree unkowns are then solved for each space curve and the developable surface.144This approach uses the same definition but only differs after the follow differentialequation has been derived.drdu ds dl- g- nds dl— g g — n gR +s +1411dty dty dtydL xdot product withy( dp’du1\ dt) dt_.; ;;— dt dt+Sdt dt +ldt dt(11.10)(H.11)Appendix H. Relating the Space Curves to the Developable SurfaceFI dux —g XLI duy —gy flyduz —g nzThe Second ApproachdudtdsdtdldtTaking the dot product with vector g:—( dudt g +s-dng +l-This simplifies to the following:ds ( ...I’duSimilarly, dot product with n,and simplifying yields:dl ( j’\du=Similarly,(11.1)(11.2)-ï (H.3)(11.4)(11.5)(H.6)(H.7)(11.8)(11.9)Appendix H. Relating the Space Curves to the Developable Surface 145Note perpendicularvector relationshipsForum = g;_dt dt — dt2 g—dt dtbut, identity:(A x B) Ctherefore, xSimphfyrng:/ —..1 .12= xg dt(;.;+xdt(CxA)B/—.-._—.—.\—s.(dp dp\dg=dtiç I== i•_ _.g+l(BXC)A(11.12)(H.13)(11.14)(11.15)(11.16)(H.17)dt(H.18)(11.19)(H.20)x g (H.21)IT(11.22)(11.23)d12(11.24)(11.25)(11.26)(&but Identity:(A x B) C(c _;therefore,—x Y).( ;du ) dtyielding:dt2— x g-—dp dp— dt dt dt2..12I P—• g —‘-• fl(dp dp d2p/ — .—.—.._.(sg +ln))du =(dr dp’\Appendix IIntersection of Developable SurfacesOnce the developable surfaces have been created there is then the requirement thatthey must intersect without gaps or spaces between them. Shown below is the vectorrepresentation of the initial conditions and their orientation relative to each other. InFigure 1.1 we see the various vectors in which their dependence to each other will beshown as follows:VFigure 1.1: Intersection of Three Developable Surfaces4Utdr1dt146Appendix I. Intersection of Developable Surfaces 147As shown in Figure 1.1 the independent parametric variable is t. The two dependentvariables are u and v. There are three developable surface control vertices with all functions of the independent variable t. The independent space curve is p(t); one dependentspace 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 r1which is relating f(u) and p (t) and the other intersection curve r2 which relates g (v) andp (t).Derivation of the equations used to calculate the intersection of the developable surfaces are as follows:= p+si: (1.1)dr1 dp dg dsit=- + si -- + -- g (1.2)dr1 (df---_xg = 0 (1.3)=(1.4)dp (df dg (dfL -— gj +s1-- Xds1 \ / (15)d2p.--16dt—dt (.)dt did2fd — — !. (17du duuju wdt dur1= f +s2ugu (1.8)Appendix I. Intersection of Developable Surfaces 148dridridu (1.9)dr1 du dg dsu) (Ho)+ S2u— + g—dtdt du du dudr df - du (dhdfu dg ds2Uu du du du— — + S2u — + (1.11)Since is perpendicular to the last term drops out (1.12)duSubstituting gives, (1.13)dr1 df du df df ‘‘ df=du— _- ( —(1.14)• — — S2uBut, (1.15)r1= f,4 + s2ugu, and (1.16)2ugu =— fu (1.17)and, r1 = p + s1tg (1.18)Substituting gives, (1.19)dffudu(dhdhd2f- (1.20)u du dudu(dh4fud2f= _.(ri_fu)) (1.21)du(dhd2f=—_.T_T.(P+s1_fU)) (1.22)(;du__________________(1.23)dgdgdu (1.24)----.i= P + (1.25)Appendix I. Intersection of Developable Surfaces 149dr2dtdr1 (dhvI\dr2 = h+siS1,j9.j = (1.39)— ; ds2—=0— dp (dh— dt idvdp Idh —•.dg+S2tdg+S2-;-.(dh —Idh -I\d xds2idtdgdtdgdtr2dr2dtdr2dtdr2 dhdt dv(1.26)(1.27)+ ( x (1.28)(1.29)(1.30)(1.31)(1.32)(1.33)(1.34)(1.35)(1.36)(1.37)(1.38)—X g,d2p — —— •iti•Yt dpdtdp dpdt dtd2h ———-•g dh-dv dv= h+sidr2 dvdv dtdv I dh dg dsi= + s + g,dv (dh dh dg dh ds2 dh=Since . = 0 the last term drops outdvSubstituting for gives,r2 = p+S2tTh (1.40)Appendix I. Intersection of Developable Surfaces 150dr2dhdv(dhvdI= . -— .- (1.41)dv(dhvd1d2h= — (r2 — h0)) (1.42)dv(dh0d’d2h= (1.43)( dgdv++ .--(1.44)— (p + S2dgdgdv- (1.45)Appendix JDerivation of Flat PlatesOnce a developable surface has been created, the shape and position of the consecutiveflat plates that make up the developable surface need to be known. Below is a simplederivation of the algorithms used to provide the information of the shape of the plates inthe x-y plane and the angle and rate of change of bending of the plates relative to eachprevious one.In order to show the relations we need to define the various vectors in the two differentspaces.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 . The vector curvatureat a point along the directrix is defined as . The vector generator at a point along the3 dimensional directrix is defined as 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 . The vector curvature, to be integrated to find the next tangent, is defined as.The vector generator ata point along the 2 dimensional directrix is defined as g2d.Since we are resolving the desired information in the 2 dimensional x-y plane, thenormal, n, is obviously n = (0,0, 1).Setting up the differential equation we intend to solve, we resolve the vector relationships into either in-plane or out-of-plane components. The following differential equation151Appendix J. Derivation of Flat Plates 152uses the existing definition of the differential equation defining the calculation of a generator and relates that as the in-plane component which defines the first term of thedifferential equation. The out-of-plane component is just resolving the out of plane curvature of the directrix and the three dimensional generator from the three dimensionalsystem.The resulting differential equation is used to find the new tangent and correspondingpoint on the x-y plane where the next point of the plate is located:2 /2dq ‘d±2 dt)dt dp= () +-g3d g2dAppendix KO.D.E. Class, Adaptive Step Size and Runge-Kutta MethodPrevious attempts at using different numerical integration methods proved that a fairamount of processing overhead takes place when making a function call to a numericalintegration 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 thanone function call is being made, like in a loop, then using an “inline” function showsconsiderable gains in speed, smaller executable files, and less functional overhead.Evaluating the Initial Value Problem (I.V.P.):y’=f(x,y) where y(x0)=yRunge-Kutta formula involves a weighted average of values of f(x,y) taken at differentpoints in the interval:Xfl X XnI1The fourth-order Runge-Kutta is written in one of the classical forms as follows:hYn+i = y, + + 2(m2 + m3) + m4)153Appendix K. O.D.E. Class, Adaptive Step Size and Runge-Kutta Method 154where m1 = f(x,y)m2 = f(x + + hmi)m3 = f(x + i-h, y, + hm2)m4 = f(x + h,y + hm3)The sum (ml+2m2+2m3+4)can be interpreted as an average slope. m1 is the slope atthe left-hand of the interval, m2 is the slope at the midpoint using the Euler formula togo from x to x + , m is a second approximation to the slope at the midpoint. Finally,m4 is the slope at Xr, + h using the Euler formula and the slope m3 to go from x tox + h.This is a very accurate formula ( halving the step size reduces the local formula errorby the factor it is not necessary to compute any partial derivatives of f.Another note should be mentioned, that if f does pj depend on y, then:m1 = f(x)1m2 = m3 = f(x + h)m4 = 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)]Simpson’s rule has an error proportional to h5 which is in agreement with error inRunge-Kutta formula [6].Appendix LNon-Linear Optimization Downhill Simplex MethodThe Downhill Simplex Method is a non-linear optimization tool for multidimensional minimization problems, ie. finding the minimum of a function of more than one independentvariable. The original paper citing this method can be found by the authors Nelder andMead [16]. A more recent overview with code samples can be found in Press [20]. Aquick overview of the method and code implemented is presented here and mentionedbriefly 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 threedimensions it is a testrahedron, not necessarily the regular tetrahedron. The algorithmis supposed to then make its own way downhill through the geometrical configuration ofan 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 initialsimplex. It then starts and takes a series of reflections, contractions and expansions. Eachone of these is assigned a variable: cv, reflection; /3, contraction; y, expansion. The firstvariable, a, is a positive constant, which is a scalar multiplication constant which mirrorsthe point through the simplex. The second variable, /3, is a constant which values liebetween 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 thecentroid 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.155Appendix L. Non-Linear Optimization Downhill Simplex Method 156Illustration 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 MethodA few trial runs were done varying the three variables, c,,-y and are listed in thetable L.1 below.*define ALPHA 0.9*define BETA 0.4*define GA11IA 1.9*define GET_PSUW for (j0;j < ndi.;j++) { for (i0,sii.0.0;i<mpts;i++)\su + pLi] [j]; psu[j]su;}Appendix L. Non-Linear Optimization Downhill Simplex Method 157a/37202 42 13 31.0 4.01.0 2.00.9 0.4 1.9Table L.1: Values of variables used for Simplex Methodvoid si.plex(double **p,double *y,iat ndia,double ftol,double (efunk) (double *),hit *nfunk, hit nan){mt i,j,k,ilo,ihi,hthi,apts—ndiatl,rnnout;mt 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_PSUMfor (;;) {iloO;ihi = y[O]>y[1J ? (inhil,O) (inhiO,1);for (i0;i<apts;i++) {if (y[i] < ylio]) ioi;if (yti) > ytihi)) {inbiihi;ihii;}else if (y[i] > y[hthi))Appendix L. Non-Linear Optimization Downhill Simplex Method 158if (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<<”\nlif( 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] [JU;ph) [j]psla[j);}y[iJ(*funk) (psusO;if(kbhitO) goto runout;}}*nfunk + ndia;GET_PSUMAppendix L. Non-Linear Optimization Downhill Simplex Method 159}}:irunout: 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 BETAtundef GAMMAAppendix MCodelistingM.1 Class Tools UsedM.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);160Appendix M. Godelisting 161vectork 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 operatorfriend vector operators(const matrixi a, coast vectort b);vectort operator*(float a);vectort operatorn(const vectort a); II elaent—by—eleaent multiplicationfriend vector operator/(const vectort v, float a);friend istreant operator>>(istreaak, vectort);friend ostresat operator<<(ostresat, coast vector&);tendifM.1.2 matrix.htinclude <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 162friend mt 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 NULL);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&);#endifM.1,3 develop.htinclnde ‘matrix .h”tinclude <fstream hpp>*ifndef_CURVE.#defInc _CURVE_class Curve{hit n;char splinetype;vector *point;static matrix BB;public:void aptype (char what) {splinetypenhat ; };char typeofsplineO{return splinetype ;};Appendix M. Codelisting 163tnt 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, mt i);double deriv_c(double t, 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 164void Draw_2d(iat genno=5O,double steop7.O,int surfao=1);surfaceO;Surface(Surface& s);Surfacet operator(Surface& s);Sendiftifudef _FILELISTS_Idefhte_FILELISTttinclude <fstreaa.hpp>class Filelisto{public:ofstreaa file;Filelisto *nezt;ofstreaat operatorO(int tO;class Filelisti{public:ifstream file;Filelisti tusit;ifstreamt operatorlj(int n);tendifM.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);Appendix M. Codelisting 165public: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);void rk(double to, doublet del_t, double ti, vectort z,vector (*f)(double t, vectort x), vectort err_scale);tendifM.1.5 vector.cpptinclude <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 166{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);Ifor(int 1=0; i<n; j++)Appendix M. Codelisting 167elaent [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;Appendix M. Codelisting 168for(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];Appendix M. Codelisting 169return 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;Appendix M. Godelisting 170element = new float * [nr];for(i0; i<nr; i++){element[i] = new float[nc];for(j0; j<nc; j++)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] = new floattnc];element[i][O] = 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++)Appendix M. Codelisting 171clement [ii [ii = a. element [i] [ii;}return *this;}vector matrix: :row(int i){vector out(nc);for(int j0; j<nc; j+)out[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 172aatriz 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] Li];}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;Appendix M. Godelisting 174Ireturn teap;Isatrix 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;Ireturn 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);IIreturn(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] [ii <0.0)teap[iJ[j) = —teap[i][jJ;IIreturn te.p;Appendix M. Codelisting 175}M.i.T solve.cppS include <fstream . hpp>Sinclude <stdlib.h>tinclude “matrix .h”coast double TINY = i.Oe-35;inline void error(char * message) {cerr < message;exit(1);};inline double abs(double & a) {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 largestdouble 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 columnsfor(i0;i<j;i++){ I/solve for the elements of the upper triangularAppendix M. Codelisting 176double sum = c.elementLi][j]; If matrix.for(int k0;k<i;k++)sun —= c . element Li] [k] *c . element I:k] U];c . element Li] [j] = sum;}double largest = 0.0;mt imaxj;for(ij;i(a.nr;i++)( I/solve for the elements of the lover triangulardouble sum = c.elementfi]Lj]; I/matrix.for(int k0;k<j ;k++)sum —= c . element Li] [k]*c .element Uk] Li];c . element 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){ f/if necessary, interchange rows.ptemp = c.elementLj]; //IOTE: this row interchange method depends onc.elementLj] = c.element[imsx]; f/the method of storing the matrix.c.elementLimax] = ptemp;ptemp = b.elementLj]; f/Interchanging rows of the RES matrix nowb.elementLi] = b.elementLimax]; I/relieves the necessity of tracking theb.elementLimax] = ptemp; f/row interchanges for later use.if(det!=NULL)*det = (*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 177for(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 = b[i] [j);mt 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++)sum —= c.element[i][k]tb[k][j];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 178for(int i0;i<a.nr;i++){ i/loop over all rows, finding the largestdouble 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 (largest 0 .0)return 0.0;scale Li) =1.0/largest;}for(int j=0;j<a.nr;j++){ //loop over all columnsfor(i0;i<j;i++){ i/solve for the elements of the upper trianguJLardouble sum = c.elementLi)Lj); // matrix.for(int k0;k<i;k++)sum—= c .elementLi) Lk]ec.elementLk) Li];c.elementLi)[j) = sum;}double largest = 0.0;mt i.marj;for(ij;i<a.nr;i++){ i/solve for the elements of the lower triangulardouble sum = c . element Li) U); //matrix.for(int k0;k<j ;k++)sum—= c . element Li) Lk) Cc. element Uk) U);c.elementUi)Uj) = sum;double tempscaleLi)eabs(sum); //keep track of the largest element.if(temp>largest){largest = temp;imaxi;}}if(imaxj){ //if necessary, interchange rows.ptemp = c.elementLi); /iIOTE: this row interchange method depends onc.elementUj) = c.elementLimax); //the method of storing the matrix.c.elementLimax) = ptp;deter(deterl .0)?—1 .0:1.0;scale Lims-x)scale Li);}Appendix M. Codelisting 179if(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 natrix.deter * 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 Iftinclude <nath.h>tifdef __ZTC__#include <fstrean.hpp>telse*include <fstrean.h>SendifSinclude “vector .h”finclude “natrix .hSinclude “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){Appendix M. Codelisting 180nasplinetypea. 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;Appendix M. Codelisting 181tt[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 == ‘2’) {if(t2 == 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 182pi[3][O] = ((point[t2+1][O])*(1.O— bb2)— (point[t2] [O])*bbl— (point[t2—1] [O])tbbO)/bb3;pi[3][ ] = ((point[t2+i][1])s(1.O— bb2)— (point[t2] [i])*bbi— (point[t2—i] [1])*bbo)/bb3;pi[3][2 = ((point[t2+1][2])*(1.O — bb2)— (point[t2] [2])*bbi— (point[t2—1] [2])*bbo)/bb3;}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 183pi[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 —t3 = 1.0;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];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] = point (t2—i] (0];p1(0] [1] = point (t2—i] [1];p1(0] (2] = point [t2—i] [2];}if(t2 = n— 2) {pi[3][O) = ((point[t2i-1]tO])s(i.O— bb2)— (point [t2] [0] )*bbi— (point(t2—i] [0])sbbO)/bb3;pi[a](i] = ((point[t2+i](i])s(i.0 — bb2)— (point[t21 [i])*bbi— (point[t2—i] (i])*bbo)/bb3;pi(3](2] = ((point(t2+i](2])s(i.0— 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 186tang[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 == n —t3 = 1.0;t2 = 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];}Appendix M. Codelisting 187if(t2 == n— 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— bb2)— (point [t2) [0))sbbi— (point[t2—1] [0])ebbo)/bb3;pi[3][1] = ((point[t2+i][i])*(1.0— bb2)— (point[t2] [1])sbbi— (point[t2—1] [1))sbbO)/bba;pi[3][2] = ((point[t2+i)[2))*(1.0— bb2)— (point[t2] [2])*bbiAppendix M. Codelisting 188— (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){pi[0] [0] ( (point [t2] [0]) *2— (point [t2+1] [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];Appendix M. Codelisting 189pi[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 == a—t3 = 1.0;t2 = 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 190else 11(1 == t2) den, = 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,1992else 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, 1992else if(i == t2 + 2) deny = cc[0][3];else deny = 0.0;}else if( splinetype == ‘1’){jf( 1 = t2— 1) deny = cc[0][0];else if(i = t2) deny = cc[O][1] +(t20.O?cc[O][0]:O.O);else if(i == t2 + 1) deny = cc[0][2] + (t2n—2?cc[0] [3] :0.0);else itCi == t2 + 2) deny = cc[0][3];else deny = 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— 1) deny = cc[0][0];else if(i == 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 == t2 + 2) deny = cc[0][3];else deny = 0.0;}return deny;}double Cunve::deniv..c(double t, mt 1){jut t2;Appendix M. Codelisüng 191double t3,deriv;atrix tt(1,4),cc(1,4);t2(int)floor(t);t3t—(double)t2;if(t == n —t3 = 1.0;t2 = — 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 = cc[0][1] +(t20.0?2.0*cc[O) [0] :0.O)+(t2u—2?—cc[O] [3] :0.0);else if(i t2 + 1) deny = cc[0J[2] + (t20.0?—cc[0)[0]:0.0)+(t2u—2?2.0*cc[0][3]:Q.O);else if(i == t2 + 2) deny = cc[0][3);else deny = 0.0;}else if( splinetype == ‘1’){if( i == t2 — 1) deny = cc[0][0J;else if(i t2) deny = cc[0][1] +(t2O.O?cc[Q][O]:O.O);else if(i == t2 + 1) deny = cc[0J[2] + (t2n—2?cc[0] [3] :0.0);else if(i == t2 + 2) deny = cc[0][3);else deny = 0.0;}else{if( i t2 — 1) den, = cc[0][0];else if(i == t2) den, = cc[0][1] +(t2O.O?2.O*cc[O][O]:O.O)+(t2—2?—cc[O][3]:O.O);else if(i == t2 + 1) deny = cc[0J[21 + (t20.0?—cc[0]f0):0.0)+(t2u—2?2.0*cc[0][3]:0.0);else if(i == t2 + 2) deny = cc[0][3];else deny = 0.0;}Appendix M. Codelisting 192return dent;}void Curve: :initbasis(char spitype, double bi, double b2){double del;if(spltype ‘1’){BB[0][0) -bi;BB[0][1) 2.0—bi;BB[0][2] = bl—2.0;BB[0][3] bi;BB[1] [0] = 2.0*bl;BB[1][1] = bl—3.0;BB[1][2] = 3.0—2.0*bl;BB[1][31 = —bi;BB[2][0] —bi;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] = 6.0*((bl)*(bl)*(bl))/del;BB[1][1] = 3.0*((b2)+2.0*((bl)*(bl)*(bl))+2.0*((bl)*(bl)))/del;BB[1][2] 3.0*((b2)+2.0*((bl)*(bl)))/del;BB[1][3] = 0.0;BB[2][0] = —6.0$((bl)*(bl)*(bl))/del;Appendix M. Codelisting 193BB[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) = 2.0*(bl)*(bl)*(bl)/del;EB[0] Ii) = 2.0t((b2)+((bl)*(bl)*(bl))+((bl)*(bl))+(bl))/del;BB[0] 12] = —2.0*((b2)+(bl)*(bl)+(bl)+1.0)/del;BBEO][3] = 2.0/del;BB11]10] = 6.0*((bl)*(bl)*(bl))/del;BB11][1] = 3.0*((b2)+2.0*((bl)*(bl)*(bl))+2.0*((bl)*(bl)))/del;BBI1)[2) = 3.0*((b2)+2.0*((bl)S(bl)))/del;BB[1][3] = 0.0;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 = te.p—>aext;return te.p—>file;}Appendix M. Codelisting 194ifstreaat Filelisti: :operatorfl (tnt 11){Filelist i ttespthis;for(; n>O; n--) te = te.p->next;return teap—>file;}M.1.9 ode.cpptiuclude “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 195{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;}IndexAerospace, 2Alignment of End Generators, 33alignment process, 36Allowment of Different Number of Control Vertices, 48Arctic Fishing Vessel, 100Areas of Application, 1AutoCAD, 96Barsky, 17Basis functions, 16C++ Classes, 85CAD Program Environment, 96CAD Program Selected, 96Change in Angle With Respect to UnitMotion of a Control Vertex, 34Class Curve, 91Class matrix, 90Class ODE, 93Class Surface, 91Class vector, 85Clement, 9Code Portability to Different Platforms,85Computer Language Chosen, 83concept of mobility, 38Conical Type Surface, 77Constraints Defining Modified Conventional Approach (Modified Nolan’sApproach) in Order to Create aNormal Directrixlan’s Approach) in Order to Create aNormal Directrix, 42continuity, 18Controlling Alignment with Mobility, 39Creation of a Normal Directrix from TwoSpace Curves, 30Demonstration Examples, 98dependent parametric variables, 73Derived Equations Yielding Intersectionsof Surfaces, 74develop.cpp, 179develop.h, 162Developable Mobius Strip, 98developable mobius strip, 39developable surface, 4directrix, 4196Index 197Downhill Simplex Method, 67Dunwoody, 10Dunwoody and Konesky, 31Equations Yielding Approximated Normal Directrix, 52explicit non parametric form, 14F117A stealth fighter, 2fanning, 49flat plate equation, 81Flat Plate Layout, 80flat plate layout, 81geometric continuity, 19in-plane curvature, 46independent parametric variable, 73Integral Chosen to Minimize, 64Interplate Angles, 80Intersection of Developable Surfaces, 73Intersection of Developable Surfaces andFlat Plate Layout, 73Kilgore’s Manual Graphical Solution, 6least squares approximation, 64Manufacturing, 2matrix.h, 161Methodologies, 5Mobius Strip, 38Modern Approach Utilizing a Single Normal Directrix, 31Naval Architecture, 1new approach, 12Nolan, 8non parametric implicit formulation, 15non-equal number of control vertices, 46non-uniform curves, 18Normal Directrix Approach, 10normal directrix control vertices, 50ode.cpp, 194ode.h, 164Offset of Normal Directrix, 47OOA- Object-Oriented Analysis, 84OOD- Object-Oriented Design, 84OOP- Object-Oriented Programming, 84Open Architecture, 96Optimization Function, 60Optimization of a Normal Directrix, 60out-of-plane curvature, 46parametric continuity, 18parametric curves, 15matrix.cpp, 169Index 198PC 386/486 Environment, 95 Uniform Non-Rational Beta-Spline, 24phantom, 78 Uniform Non-Rational Tension CatmullPhantom Surfaces, 78 Rom Spline, 21phantom surfaces, 77 Utilizing Modified Conventional ApproachPNA, 6 to Approximate a Single Directrixrate-of-rotation of a generator with reDirectrix, 50spect to motion of the control vertices along the normal directrix vector.cpp, 165ntrol vertices along the normal direc- vector.h, 160trix, 34rational continuity, 20ruled surface, 3-Simple Conical Developable, 100Simplex Parameters Used, 64small in and out-of-plane curvature, 41solve.cpp, 175Splines, 14Textiles, 3The Downhill Simplex Method, 68threshold, 42tolerance, 42Types of Splines Chosen, 14UBC Series Fishing Vessel, 100uniform curves, 18
- 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 Issued | 1993 |
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 |
FileFormat | application/pdf |
Language | eng |
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. |
IsShownAt | 10.14288/1.0080847 |
URI | http://hdl.handle.net/2429/4930 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1994-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1994-0078.pdf [ 3.56MB ]
- Metadata
- JSON: 831-1.0080847.json
- JSON-LD: 831-1.0080847-ld.json
- RDF/XML (Pretty): 831-1.0080847-rdf.xml
- RDF/JSON: 831-1.0080847-rdf.json
- Turtle: 831-1.0080847-turtle.txt
- N-Triples: 831-1.0080847-rdf-ntriples.txt
- Original Record: 831-1.0080847-source.json
- Full Text
- 831-1.0080847-fulltext.txt
- Citation
- 831-1.0080847.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
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