Fast Contact Evolution for Piecewise Smooth Surfaces by Paul G. Kry B.Math., University of Waterloo, 1997 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Department of Computer Science) We accept this thesis as conforming to the required standard The University of British Columbia August 2000 © Paul G. Kry, 2000 In presenting this thesis/essay in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying for this thesis for scholarly purposes may be granted by the Head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. 2PCO io m Date Department of Computer Science The University of British Columbia 2366 Main mall Vancouver, BC Canada V6T 1Z4 Abstract Dynamics simulation of smooth bodies in contact is a critical problem in physically based animation and interactive virtual environments. We describe a technique which uses reduced coordinates to evolve a single continuous contact be-tween Loop subdivision surfaces. The incorporation of both slip and no-slip friction into our algorithm is straightforward. The dynamics equations, though slightly more complex due to the reduced coordinate formulation, can be integrated eas-ily using explicit integrators without the need for constraint stabilization. The use of reduced coordinates also confines integration errors to lie within the constraint manifold which is preferable for visualization. Our algorithm is suitable for piecewise parametric or parameterizable sur-faces with polygonal domain boundaries. Because a contact will not always remain in the same patch, we demonstrate how a contact can be evolved across patch bound-aries. We also address the issue of non-regular parameterizations occurring in Loop subdivision surfaces through surface replacement with n sided S-patch surfaces. Three simulations show our results. We partially verify our technique first with a frictionless system and then with a blob sliding and rolling inside a bowl. Our third simulation shows that our formulation correctly predicts the spin reversal of a rattleback top. We also present timings of the various components of the algorithm. ii Contents Abstract ii Contents iii List of Figures vi Acknowledgements viii 1 Introduction 1 1.1 Chapter Summary 3 2 Related Work 6 2.1 Dynamics Simulation ." 6 2.2 Contact Kinematics 9 2.3 Smooth Surface Contact, Collision and Distance 10 3 Additional Background 12 3.1 State Representation 13 3.2 Constraints 13 3.3 Solving Constrained Systems 14 3.4 Reduced Coordinates for Constrained Systems 15 4 Formulation 19 4.1 Preliminaries 19 iii 4.2 Contact Kinematics 21 4.2.1 Contact Coordinates 21 4.2.2 Multiple Contacts 22 4.2.3 Kinematics Equations 24 4.3 Constrained Dynamics 27 4.4 Friction 29 4.4.1 Dynamic Friction 29 4.4.2 No-Slip Friction 32 4.4.3 No-Slip and No-Spin Friction 34 4.4.4 Changes in Friction . . • 35 4.5 Traversing Patch Boundaries 35 4.5.1 Time of Crossing 38 4.5.2 Contact Location on Adjacent Patch 39 4.5.3 Computing New Contact Coordinate Velocities 40 4.5.4 Corner Transitions 41 4.5.5 Avoiding Traversal Problems 42 4.6 Algorithm Summary 43 5 Surface Representation 45 5.1 Subdivision Surfaces 45 5.1.1 Parametric Evaluation of Subdivision Surfaces . 48 5.1.2 Non-regular Parameterization 49 5.2 Surface Replacement 51 5.2.1 S-patches 52 5.2.2 Matching Continuity 56 5.3 More on Computing Derivatives and Blossoms 62 5.3.1 Derivatives of S-patches 64 5.4 Implementation Description 64 iv 6 Simulation Results and Evaluation 67 6.1 Timings 67 6.2 Bowl on Flat Frictionless Surface 69 6.3 Blobs in Bowls 72 6.4 Rattleback Simulation 73 7 Conclusions 77 7.1 Future Work 78 Bibliography 79 Appendix A Contact Kinematics Derivation 84 Appendix B Example Files 87 B . l Scene Description File 87 B.2 Subdivision Surface File 89 v 5 List of Figures 1.1 Example of dynamic simulation 1 3.1 Planar systems with reduced degrees of freedom 15 3.2 Four possible configurations of a five bar mechanism 17 4.1 Contact Coordinates 23 4.2 Change of reference frames 25 4.3 Two bowls showing patch boundaries 36 4.4 Domain boundaries and functions involved in a transition 37 4.5 Location of intersection 40 4.6 Angle correction for ip at a patch boundary 41 4.7 Transition through a corner 42 5.1 Loop subdivision rules 46 5.2 Three levels of subdivision of an octahedron 47 5.3 A planar valence 3 Patch 50 5.4 A planar valence 10 patch 50 5.5 An S-patch fit to an extra-ordinary vertex of a bowl 52 5.6 S-patch maps 53 5.7 Adjacent domains 57 5.8 Middle refinements used to create an S-patch 58 5.9 Non uniquely determined control points 59 vi 5.10 Tangent plane at a corner of an S-patch 60 5.11 Code produced with automatic differentiation techniques 65 6.1 Bowl on a flat frictionless surface 70 6.2 Graph showing the position of the center of mass over time 71 6.3 A lumpy marble rolling on a plate 72 6.4 Two different rattleback tops with asymmetric shape 73 6.5 Our rattleback model, top and front views 74 6.6 Tap start of a rattleback 75 6.7 Spin start of a rattleback 75 6.8 Sequence showing a second spin reversal 76 vii Acknowledgements Many thanks go to my supervisor, Dinesh K . Pai, as this thesis would not have been possible without his guidance. I would also like to thank Uri Ascher for his invaluable comments and suggestions. Lastly, to all the others who are too numerous to name here, I am very grateful. This thesis was supported in part by scholarships from N S E R C and BC-ASI . PAUL G . KRY The University of British Columbia August 2000 viii Chapter 1 Introduction Within the graphics community, interest in physical simulation of rigid bodies prob-ably reached a peak sometime in the last ten years. Although some might consider it to be a solved problem, we feel that there are still important challenges to be met. One of the difficult aspects of simulating rigid bodies is the case where there is contact. In this case the motion of the bodies must be constrained so that they do not in-terpenetrate. The combination of smooth surfaces, friction and multi-ple contacts can make this a very dif-ficult problem. This thesis is about the dv- r- < 1 I c J I rigure 1.1: Example ol dynamic simula-namics of piecewise smooth surfaces tion using our algorithm. The objects are bounded by smooth Loop subdivision sur-m single continuous contact. We c ° laces. present a method of computing how a system of two rigid bodies evolves when their smooth boundaries are in contact at a single point. We call this problem the contact evolution problem. Figure 1 shows 1 an example of a smooth globe-like object rolling and sliding inside a bowl, simulated with our system. We formulate the problem as an ordinary differential equation in the coordi-nates of the contact patches. By parameterizing the system in reduced coordinates, violation of the contact constraint is prevented as we integrate. Any errors that accumulate during integration instead lie entirely within the constraint manifold. The incorporation of dynamic Coulomb friction into our method is straightforward. Additionally, we can easily reformulate the equations to include no-slip and no-spin friction. By monitoring the constraint force we can accommodate changes in friction and object separation. This is an attractive approach to the problem because it avoids the draw-backs of traditional methods which spend most of their time in collision detection and constraint computation. It is important to recognize that although collision detection is still necessary to identify new contacts, it can be performed at a much less frequent rate because it is not necessary for our single continuous contact evo-lution technique. Other techniques rely on collision detection or minimum distance computation at each time step to track the local contact points or for constraint stabilization. Simulating realistic models constructed with subdivision surfaces or free-form parametric surfaces such as N U R B S is important because of the popularity of these surface representations for modelling and design. These types of piecewise parametric and parameterizable surfaces are ideal for our single contact evolution technique. With respect to piecewise parametric surfaces, an important contribution of this thesis is a technique to evolve a contact across a patch boundary. The surface description of any object will almost always consist of a collection of several patches. Part of the popularity of subdivision surfaces is attributable to the relative ease with which they can be edited and implemented. Unfortunately, this facility does not carry over to contact evolution. Almost every object defined by a subdivi-2 sion surface will have multiple patches that are non-regular (in the natural param-eterization). In this thesis we focus on the use of Loop subdivision surfaces. For contact evolution to work we must avoid evolving the contact near non-regularities. We achieve this by replacing portions of the Loop subdivision surface with n sided S-patches which are more complicated but better behaved. The replacement sur-faces we construct maintain tangent plane continuity at their boundary with the original subdivision surface. As a result their shape often matches quite closely that of the non-regular surface. 1.1 Chap te r Summary We provide an account of related work in Chapter 2. The related work generally falls into the three categories: dynamics simulation, contact kinematics, and collision detection. Chapter 3 provides additional background and motivation. We start with a discussion on the representation of state. Sections 3.2 and 3.3 introduce the concept of a constraint and briefly discuss solving constrained systems. This leads into Section 3.4 where we motivate the use of reduced dimension parameterizations for Constrained Systems. Chapter 4 consists of a detailed derivation of our algorithm. After some pre-liminaries, contact kinematics equations are derived in Section 4.2 and then incor-porated as a constraint into the Newton-Euler equations of motion in Section 4.3. With the foundation of our technique explained in the first part of this chapter, we continue by demonstrating how to add Coulomb friction to the system in Sec-tion 4.4. We also present some ideas for implementing no-slip and no-spin friction in our framework as well as ideas on switching between different types of friction. Section 4.5 contains an explanation of how a contact traverses the boundary of a patch into an adjacent patch. Traversing boundaries is important as most models consist of a collection of parametric surfaces. Lastly, we conclude this chapter with 3 a summary of the algorithm. In Chapter 5 we examine the problems involved with using Loop subdivision surfaces with our algorithm. After a brief explanation of Loop subdivision surfaces, we discuss the standard method for evaluating Loop surface patches and the impli-cations of a non-regular parameterization on our contact evolution algorithm. In Section 5.2 we introduce the idea of replacing parts of the surface with S-patches. This section provides a description of S-patches and explains how to set S-patch control points to maintain surface continuity via blossomed Loop patch evaluations. In Section 5.3 we provide additional detail on how we compute derivatives and blos-soms. Since derivatives of S-patch functions are quite complicated we write this code with inspiration from automatic differentiation techniques. We wrap up this chapter with a brief overview of our Java implementation and its features. Chapter 6 shows the results of our algorithm. As part of the evaluation of our algorithm, we measure how long different portions of the calculation take in Section 6.1. The timings show that a large proportion of the computation time is spent in surface evaluations and surface derivative evaluations which suggests that our algorithm is best suited for surfaces which are easy to evaluate. In Section 6.2 we partially validate our implementation with a Loop subdivision surface on a flat and frictionless Bezier patch. Section 6.3 demonstrates two Loop subdivision surfaces in contact with dynamic friction. Section 6.4 finishes this chapter by showing that our implementation correctly predicts the spin-reversing behaviour of a rattleback top. Chapter 7 provides our conclusions and summarizes the impact of our algo-rithm. We also describe possible future work. Appendix A describes an alternate and more useful derivation of the contact kinematics equations. This derivation starts by only considering half of the contact and takes the least complicated expressions for describing this half of the equations. The other half of the equations are very similar and only require transformation into the same coordinates as the first half. 4 In Appendix B we provide some sample files and a quick description of other data files needed for our implementation. A rattleback scene description file is provided in Section B . l and a object file for a Loop subdivision surface is provided in Section B.2. 5 Chapter 2 Related W o r k In recent years there has been a good amount of work concerning simulation of contact. Being such a broad topic, the emphasis falls into several different areas. Papers coming from the robotics community tend to focus on the constraints and numerical methods while papers coming from the graphics community tend to focus on collision detection. We break the following short literature survey into three sections to help shed light on which areas our algorithm falls. 2.1 Dynamics Simulation Multibody dynamics has been a bountiful area of research for many people over sev-eral decades. Numerous books have been written on this topic such as [15, 43, 38] to give an example of just a few. It is important for us to mention this work because our research has been largely influenced by these multibody algorithms. In particular, our algorithm is inspired in part by the analysis of rigid body methods by Ascher, Pai and Cloutier appearing in [4]. Note that reduced coordinate formulations such as the one we use in this thesis are well studied (for example, see [40] by Rabier and Rheinboldt). In the rest of this section we provide an account of recent algorithms and implementations relating to dynamics simulation. Newton, [11], as described in Cremer's PhD thesis is a system architecture for 6 general purpose physical system simulation. His modular system incorporates colli-sion detection with a variety of dynamics formulations consisting of open and closed loop kinematic chains with a variety of inter-object constraints. The simulation is based on events such as collisions and changes in contact and friction. Lubich et al. present M E X X in [28]. M E X X is software for integration of constrained mechanical multi body systems. Its focus is on the numerical method. It is based on a method which is explicit in the differential equations and linearly implicit in the nonlinear constraints. It has the ability to detect events such as im-pacts, however, unlike Cremer's system it is the user's responsibility to reformulate their equations and constraints to deal with such events. Sauer and Schomer in [46] describe G A L I L E O , a system for rigid body sim-ulation with emphasis on the simulation of unilateral constraints for virtual reality applications. They formulate the problem of n polyhedral or curved objects (their examples are all portions of spheres or cylinders) with K contacts with friction in a single linear complementarity problem (they use a polyhedral approximation of a friction cone). They demonstrate their friction modeling though the simulation of a tippe-top. Stewart and Trinkle in [53] present an interesting time-stepping method for rigid body simulation. Being based on impulse-momentum equations it does not need to explicitly resolve impulsive forces. It also does not require explicit collision detection and can handle simultaneous impacts. Since they use an impulse based method rather than a force based method they can set up their problem as a linear complementarity problem. Because of their formulation, it is not necessary to de-termine a time of impact (instead they maintain a set of constraints which are in danger of being violated). They use a polyhedral approximation of a friction cone. Hahn in [19] is an example of earlier work on rigid body simulation with polyhedral models. His work focuses on integrating physically based simulation with objects which are animated by hand. Although his system uses geometric 7 constraints as well as more general non-holonomic constraints, Hahn achieves non-interpenetration by modelling contact as a series of frequently occurring collisions. BarafF's [5] is another example of some earlier work. The focus of this paper appears to be the formulation of resting contact constraints. He uses a linear com-plementarity problem to determine when contacts are breaking. This paper does not address friction and the the examples given are planar consisting of polygons. BarafF, in [6], also looks at rigid body simulation with curved surfaces. A l -though the paper focuses on implicit surfaces a method of using parametric surfaces is described in one of the appendices. Central to the method is the tracking of all pairs of extremal points, that is the points where maximum penetration occurs between the surfaces. To solve multiple contacts they formulate the unilateral con-straints as a linear complementarity problem (LCP) . Here Lemke's algorithm for solving linear complementarity problems is preferred to using a heuristic solution to the positive semi definite (PSD) linear programming problem. Baraff points out, however, that the problem is non-PSD in the presence of friction and suggests a heuristic method would be useful in this case (which he presents in [8]). The paper also discusses a fast method for detecting collisions using cached "witness" separat-ing planes. He also discuses a method for quickly and accurately determining the collision time as the algorithm requires this for backing up to the times at which collisions occur. In [8], Baraff shows how to compute contact forces with friction for rigid bodies. His approach in this paper is to generalize the problem as little as possible (he does not formulate an L C P as before). This has the advantage of avoiding code for solving general optimization problems which can be overly complicated for the friction problem. The result is a much faster and simpler algorithm. BarafF admits that although they can not prove that the algorithm terminates they have not experienced this problem from any of their tests. They implemented an interactive planar algorithm, and an off line system For the three dimensional version oF the 8 algorithm. Anitescu, Cremer and Potra's [1] provides an example of more recent work in rigid body dynamics. This work focuses on the formulation of contact problems with friction as a linear complementarity problem. The friction cone is approximated by a polyhedron. They consider smooth shaped bodies and derive the contact kinematics equations for surface-surface contacts, curve-surface contacts, and curve-curve con-tacts. They also investigate conforming planar contacts with arbitrary boundaries. They give the example of a cylindrical object on a table with external forces to demonstrate the transition from the conforming planar contact to a curve-surface contact. Although the Newton-Euler formulation they use is in global coordinates, instead of using an optimization algorithm for determining the new contact points they use the contact kinematics equations. They do not address the problem of changing patches as the contact moves since they assume that all calculations are done within a single patch on each body. Mirtich and Canny take an interesting approach to rigid body dynamics simulation in [29]. A l l interactions between the bodies, from colliding contact to all varieties of continuous contact are treated as collisions. Their method has the advantage of simplicity as a wide variety of systems can be simulated, including those that are difficult to simulate with constraint based methods. The disadvan-tage however is that it relies greatly on the collision detection algorithm to provide information about all collisions with accuracy. One of the many examples they use to test the algorithm is a rattleback top. 2.2 Contact K inema t i c s Smooth surfaces in contact have been extensively studied in both graphics and robotics. We've already mentioned Anitescu's work with the equations of contact, but there are several others which should be mentioned. Many of the earlier papers derive general contact kinematic equations but then demonstrate them with simple 9 surfaces such as spheres. Montana derived contact kinematics equations for orthogonal nets in [30] and showed how they could be used for fine grip adjustment and other robotics ap-plications. The beauty of Montana's derivation is that if the contacting surfaces are described by orthogonal networks then the kinematics relationship can be described with matrices from differential geometry. See also [33, 10, 25, 22] for alternate derivations and applications. Goyal derives various constraints in [17] for the purpose of smooth surface simulation. The equations in this paper were derived to allow smooth surface inter-action for Cremer's Newton dynamic simulator. The provided example interactions are limited to relatively simple objects such as a cylinder rolling on a block without slipping and a torus rolling on a plane. An alternate version of the contact kinematics equations is derived by Nelson Johnson and Cohen in [35]. They use a separation distance as one of their contact coordinates allowing them to easily measure penetration distance during contact evolution. In this paper the penetration distance is used to compute a contact force in a haptic simulation of N U R B S surfaces in contact. 2.3 Smooth Surface Contac t , Co l l i s ion and Dis tance Hegron computes rolling motion of a convex object on a parametric surface using a prediction-correction schema in [20]. He demonstrates a sphere rolling without sliding on a collection of bi-cubic Bezier patches. The sphere is placed in contact with the surface by lowering it vertically using an iterative method. They assume that the parametric surface patch grid has C 1 continuity, and that the rolling rigid object is convex. A single point of contact is maintained during the animation. Snyder examined the problem of placing smooth surfaces in contact in [49]. His algorithm allows the user to place smooth objects in stable non-interpenetrating configurations given a set of external forces. Shapes can be non-convex and place-10 merit of complicated shapes can occur at interactive rates. The key is that the objects do not behave as physics would dictate. The example of placing a spoon in a bowl shows how this is useful for object placement. Trying to place a spoon in a bowl using a physically based animation algorithm may be difficult as the spoon may bounce out before the system comes to rest. In an earlier paper, [50], Snyder and others at the California Institute of Technology investigated animating smooth surfaces. Their focus was on collision detection. Surfaces can be parametric or implicit and the shapes can vary over time. One of the key features of their solution is the generation of a set of uniformly distributed points on the interface of a conforming contact between smooth surfaces. •Their method although slow is very powerful. The paper focuses solely on collision detection and does not address how to handle the collision response and contact evolution problems. Johnson and Cohen, in [23], present a framework for minimum distance com-putations between parametric surfaces. Their approach is based on a lower-upper bound tree which allows them to quickly prune bounding volumes which are not involved in the minimum distance. They also show that their solution converges relatively quickly which is useful for time-critical applications. In Nelson and Cohen's [34], constraints in Cartesian coordinates are used for interaction with smooth surfaced mechanical assemblies. They report that their approach has a compact implementation due to the re-use of constraints. They demonstrate their system for several small assemblies including a Stewart platform. 11 Chapter 3 A d d i t i o n a l B a c k g r o u n d This chapter provides additional background on rigid body dynamics and simulation issues in order to help motivate our work. We will show why we feel a solution to the single continuous contact problem is necessary. Additionally, the discussion should help give insight into the difficulty inherent in this simulation problem and related ones. An important focus of our discussion will be on the choices for parameterizing the state of a system. This has a large impact on formulating dynamics equations. We will also discuss different types of constraints, and methods for solving the constrained equations of motion. This extra motivation is helpful because our algorithm is not easily compa-rable to other existing techniques. This is because in focusing on single contact we avoid issues which are often emphasized in other algorithms. For example, friction is a well studied problem and is one that becomes very difficult when there are mul-tiple contacts between two bodies (Baraff has shown that solving multiple contacts with friction is N P hard in [7]). 12 3.1 State Representat ion Rigid body simulations require a method for parameterizing the possible positions and motions of bodies. An object free to move in space has six degrees of free-dom. Three of these degrees of freedom are translational, while the other three are rotational. Euler angles, for example, can be used to parameterize the rotation, although not without singularity problems. A higher dimensional parameterization is often used as it can be more convenient. The rotational component may be represented by a 3 x 3 rotation matrix giving us a parameter space of dimension 12. But only matrices which are orthogonal and have determinant equal to 1 are rotation matrices so we must be careful that our state vector preserves this property, despite any changes we make. Alternatively unit length quaternions can be used to represent a rotation. This gives a total dimension of 7. Again, we must be careful to preserve the unit length of a quaternion, which happens to be a fair bit easier. If we make a change to the state which violates the condition by a small amount it is easy to project the state back onto the space of valid states as we can simply normalize the quaternion. See [55] for more on rigid body simulation and the issues of representing the state of the system. 3.2 Const ra in ts A constraint is an expression which describes a subset of states satisfying a certain property. It is often written in the form f(x) = 0 or f(x) > 0, where x is the state of the system. When the expression is an equality the constraint is bilateral, such as for a revolute joint. If the expression is an inequality then the constraint is unilateral, such as a non-interpenetration constraint. Holonomic constraints are essentially position constraints, or they are con-13 straints that can be integrated to position constraints. Examples of this are revolute joints and some kinds of non-interpenetration constraints. A non-holonomic con-straint, however, is one that is expressed as a velocity or acceleration constraint, and can not be integrated to get a position constraint. An example of this is a ball rolling without slipping on a table. The ball can achieve any position on the table with any rotation by some sequence of rolling and spinning. For more on constraints in rigid body dynamics see [24]. 3.3 Solv ing Cons t ra ined Systems Combining the Newton-Euler equations of motion with an algebraic constraint gives a differential algebraic equation (DAE) . Such equations are not always easy to solve (see [2]). The most common approach is to differentiate the constraints so that they are stated in terms of accelerations. The differentiated constraints provide a system which can be solved as an ordinary differential equation. We must, however, take care to ensure the solution lies on the constraint manifold. This can be done with Baumgarte stabilization [9]. This stabilization technique acts like a spring damper system connecting our solution to the constraint manifold. As this type of stabilization has some undesirable properties (see [3]) there may be a desire to use more expensive and more mathematically sound techniques. Sometimes it is easier to give a position constraint as a velocity constraint. This is the case for parametric surfaces in contact. The position constraint is essen-tially a distance function which lacks a simple closed form. The velocity constraint, on the other hand, is much easier to write using the normals of the contacting points (or the extremal points in the case of interpenetration). 14 X Figure 3.1: A robot arm with one degree of freedom is depicted on the left while a rod in contact with the plane having two degrees of freedom while in contact is depicted on the right. These are bilateral and unilateral constraints, respectively, for which we can express the configuration of the system more concisely with reduced coordinates. 3.4 Reduced Coordinates for Cons t ra ined Systems It can be advantageous to choose a set of reduced dimension coordinates for repre-senting the state of a constrained systems. A robot arm with one degree of freedom can be described by its joint angle rather than parameterizing its position as if it were a free body. Although the reduced coordinate dynamics equations are generally more expensive to evaluate, they do not require stabilization when they are inte-grated if the reduced coordinates exactly describe all possible configurations of the system. The robot arm with its bilateral constraint is an example of this situation. We might imagine a rod on a plane in a two dimensional environment as another example but with a unilateral constraint. Figure 3.1 shows a conceptual drawing of both these cases. For the robot arm case we only need to keep track of the angle. In the case of the rod, while it is in continuous contact with the surface we only need to keep track of the angle and the x position of the contact. It is useful to consider the full implications of reduced coordinates for para-metric surface non-interpenetration constraints. A non-interpenetration constraint is unilateral. It is only active when there 15 is contact and it can be complicated as it depends on the contacting geometry. Multiple contacts mean multiple active constraints. If we use higher dimensional coordinates then although we do not need to change our state representation when the constraints become active, we do need to stabilize these constraints to ensure that our solution remains on the constraint manifold. The alternative of reduced coordinates unfortunately does not extend well to multiple contacts. The constraint (or the constraint gradient in the case of some methods of solving this problem) may not be easy to compute for certain configurations of the given geometry. There will be a best way to handle each of these special contact configurations. For example, a screw in a threaded hole would best be treated as a system with one degree of freedom. A peg in a hole would be best treated as a two degree of freedom system. A cylinder on a table has three degrees of freedom until it is tipped up on edge in which case it has five. The problem of transitions between different formulations, however, is not trivial even for this relatively simple example (see[l]). Trying to use reduced coordinates for systems with a reduced number of degrees of freedom is not always possible. The system may have multiple states for a given coordinate location, or singularities where there are reduced degrees of freedom. For example consider a five bar mechanism. This mechanism is a linkage of four bars (the fifth bar is actually the base which connects the two ends of the linkage). We might try to use the planar location of the end effector to describe the state of the system. An arbitrary position can be achieved with the left "elbow" bent in or out and the right elbow can be in or out. This can be seen in Figure 3.2. There are other positions which can not be achieved at all due to the limited length of the arm. Additionally, when one of the arms is straight there is one less degree of freedom at the end effector. We might consider using two angles to parameterize the system, but we do not really avoid any of our problems. There are still configurations which can not be reached and there are multiple possible 16 Figure 3.2: Four possible configurations of a five bar mechanism created by having the left and right "elbows" in or out. The planar position of the dark circle is the configuration space. It is inappropriate to use a two dimensional state representation for this system even although it really only has two degrees of freedom. end effector positions for many settings of the two angles. It appears that a two dimensional state representation for the system is inappropriate. This seems to be fairly strong motivation not to use reduced coordinates, at least for this scenario. Let us go back to the smooth surface contact constraint as this is our primary interest. If there is only a single point of contact between two bodies, then the con-tact can be parameterized as a five degree of freedom joint (this is explained fully in Section 4.2.1). In this case the reduced coordinates uniquely describe a configura-tion. The additional advantage of these coordinates is that their dimension matches the number of degrees of freedom exactly, and the entire space parameterizes valid configurations. As a result, we can use these coordinates without having to worry about stabilization. Note that we will still accumulate errors in integrating the sys-tem but these errors will be inside the constraint manifold. These types of errors are acceptable because they are far less noticeable than interpenetration or separation errors. Global changes such as new constraints becoming active can be detected at a slower rate than the rate at which we evolve contacts. This is possibly the most important feature of our technique as collision detection is often the most computationally expensive component of physically based simulations. Collision detection can be minimized or even omitted if we know that the object shapes do not allow multiple contacts. 17 Additionally, single contact can be useful in a variety of situations, such as haptics, [35], and certain manipulation examples, [22]. Fast contact evolution is interesting for real time interaction with smooth objects where contact evolution information is used for sound synthesis, [13, 14]. 18 Chapter 4 Formulation 4.1 Pre l iminar ies In this thesis we use a notation for spatial dynamics which is similar to [15, 44, 21, 30, 32]. The motion of a rigid body i = 1 , n is described using a reference frame labeled i attached to the body. Unless otherwise noted, we will let this body ref-erence frame be located at the center of mass and with its axes aligned with the principle axes of rotation. The homogeneous coordinates of frame i with respect to another frame j is given by the 4 x 4 matrix --E. We will always use leading subscripts and superscripts to indicate frames. The homogeneous coordinates of a three dimensional vector x in frame i are denoted 'x. The homogeneous coordinates of this vector in frame j are given by left multiplying by jE. The spatial velocity <f>(j, i) describes the relative motion of frame i with respect to frame j. In (non-homogeneous) coordinates of frame i it is given by the 6 x 1 matrix *<f>(j,i) = ('u;T, V r ) T , where u is the angular velocity and v is the linear velocity of the point at the origin of frame i. Spatial forces, called wrenches, are represented as / = (/J, fJ)T where fr is the (rotational) torque and ft is the (translational) force. 19 0(j> = (4.1) It is useful to note that the spatial velocity vector parameterizes a twist that can be written in homogeneous coordinates as the 4 x 4 matrix, M v 0 0, Here, [u>] denotes the skew symmetric 3 x 3 matrix equivalent to the cross product u x , i.e., 0 0 -ujy -UJ' 0 J (4.2) and O denotes the linear expansion operation which converts the 6 x 1 matrix into a 4 x 4 matrix representation. We will let V denote the linear operation which performs extraction of the 6 parameters, thus, V M v o o, (4.3) If V is applied to a matrix whose upper 3 x 3 matrix is not skew symmetric then we can assume that V acts as a projection onto so3. The time derivative of the change of coordinates from frame i to frame j is a twist when written in coordinates of frame i, i.e., V (J , 0 = V ( ; E £ ) (4.4) Spatial velocities and spatial wrenches transform according to the Adjoint trans-formation ;Ad . Spatial velocities being contravariant quantities transform by left multiplying, J<p = -jAd l<f>. Because we write spatial wrenches as column vectors, these covariant quantities are transformed by left multiplying with the inverse trans-pose, 3f = ] A d T ' / • The 6 x 6 Adjoint matrix is defined as, ^Ad = I 6 ° | , where ]E >]e e, 20 Here 0 is a 3 X 3 rotation matrix and p is a 3 X 1 displacement. Lastly, we define the spatial cross product of (f> = ^ uT, v operator with coordinate matrix as the linear M o [v] M 4.2 Contact K inema t i c s We consider rigid bodies whose boundaries are defined by a collection of parametric surface patches. Although the surface patches will typically be polynomial or ra-tional polynomial functions this does not necessarily need to be the case. We only require the ability to evaluate the function and its derivatives at a given point. A single point of contact between two bodies will involve one patch from each body. We can ignore, for now, the fact that adjacent patches will overlap along their boundary curves as it is the location of the contact that we want to describe and if it lies on a boundary then we can choose either of the patches sharing the boundary to describe the contact location. 4.2.1 C o n t a c t C o o r d i n a t e s Suppose body 1 and body 2 are in contact. Let the shape of the contacting patches be described by the functions J c : (s,t) —>• R 3 and 2 d : (u,v) —> E 3 . Note that the functions have a preceding superscript denoting that they define the surface shape in the coordinates of their respective body's reference frame. We drop this superscript when the coordinate frame is clear from context. We use the notation c s = f The orthonormal contact frame with coor-dinates (s, t) can be defined in frame i , the body frame, as the coordinate frame with origin at c(s, t) and coordinate axes x = y = (C 3 • C S ) C f - (c s • C t)c s z = (4.5) | | ( c , - c t ) c t - ( C S - C ( ) C S 21 Coordinates are in frame 1 (superscripts were omitted). The vectors x, y and z always form orthonormal axes, and will always exist if we require our parametric surface to be regular (that is c ) S (s, t) ^ 0 and c i t (s , t) ^ 0 for all (s, t) in the domain). Also, y is orthogonal to z since it is a linear combination of c j S and c^ . It is also orthogonal to x since its inner product with ctS is equal to zero. Thus the homogeneous coordinates of this contact frame with respect to frame 1 are given by 1 1 1 x x y z c (4.6) 2 c C — (4.7) 0 0 0 1 The matrix in equation 4.6 is a function of s and t. The contact frame on body 2 is defined in the same way yielding the matrix | C E which is a function of u and v in the domain of the patch' d(u, v). Note that the c in the frame label 2c stands for contact as opposed to representing patch c of body 1. Because the contact frames are orthonormal and occur at the same location in space, they can be easily aligned by a rotation matrix (see Figure 4.1). / cost/? — sin V7 0 0^ — sin tp — cos tp 0 0 d e f 0 0 - 1 0 \ 0 0 0 1/ The matrix g£E is idempotent and is its upper 3 x 3 rotation matrix. We can now describe any contact configuration between patches c and d of bodies 1 and 2 by the 2-dimensional location of the contact in the domain of each patch along with the angle of rotation tp. Assembled in a column vector, we call these 5 independent variables the contact coordinates and we denote it q; i.e., q = (s t u v v) • (4-8) 4.2.2 Mult iple Contacts When there is more than one point of contact between two bodies the situation becomes more complex. When there are two contacts the number of degrees of 22 Figure 4.1: Contact Coordinates freedom is reduced in general to four. We run into problems, however, if we try to represent our two contact configurations as a four dimensional generalized joint. Suppose we use the parametric location of one of the contacts on both surfaces as our contact coordinates. Valid unique configurations in the state space are those where there is only one rotation of the objects about each other at the specified point of contact for which the objects establish a second point of contact. If there are additional rotations which yield a second point of contact then we have the problem of not knowing which configuration we have from the contact coordinates alone. For some special cases of contact there can. be in fact five degrees of freedom. An example of this is a sphere wedged between two flat parallel planes. In this case all rotations are possible giving us the extra degree of freedom. Although this case can be treated as a bilateral constraint (the sphere can not move in the direction normal to either contact), there are most certainly scenarios which can not be simplified so easily. Other attempts at a general parameterization of multi-contact configurations will have similar problems. Parameterizing multiple contact configurations is a difficult problem we do not propose to solve. 23 This all appears to be a good argument for avoiding a reduced coordinate for-mulation and sticking with simpler more reusable constraints that need stabilization. This avoids the problem of choosing and switching between different parameteriza-tions. Nevertheless, we still recognize the single contact problem as interesting, important, and useful. If the objects spend most of the time in a given contact state then parameterizing the local contact manifold is useful. 4.2.3 Kinematics Equations The contact kinematics equations relate the relative motion between two smooth contacting bodies to a change in the contact coordinates. Because the contact coordinates are of reduced dimension this system of equations will be a linear transformation (a function of the contact coordinates) from the 5 dimensional space of contact coordinate velocities into a 5 dimensional subspace of the 6 dimensional space of spatial velocities. Many different derivations of these equations with minor differences can be found in [30, 10, 1, 35]. For example, [30] showed that when the parametric surfaces are orthogonal networks these equations can be written using the fundamental forms of curvature of the two surfaces. For another example, if the contact coordinates are extended to include a distance along the z direction between the coordinate frames then we get a bijective transformation relating all spatial velocities to the changes in these extended contact coordinates. A derivation of this can be found in [35] where it is used for haptics. Having a separation distance is useful in haptics since a negative distance measures penetration which can be used to compute a reaction force based on a penalty method. In our simulation, however, because we use the contact kinematics equations as a constraint we do not include distance as one of the degrees of freedom. We will 24 1 E 2 2c 1c 2c E ic E 1 E Figure 4.2: Change of reference frames give our own derivation here, suitable for any type of parametric surface. A more detailed derivation can be found in Appendix A. As illustrated in Figure 4.2 we compute a matrix representing a change in coordinates from frame 1 to frame 2. This composition is written as ?E = LE £E \CE. (4.9) The three matrices on the right hand side are functions of the contact coor-dinates. They are defined by Equations 4.6 and 4.7. Recall that the relative spatial velocity can be defined by 1) = V \E ? E. (4.10) We take the time derivative by using the chain rule, thus 5 V ( M = v £ {E f E l W f l . (4.11) i=i 25 Since the contact coordinates can change independently of each other, each 4 x 4 matrix in the sum will be a twist and can have the extraction operator applied directly, hence 5 1<t>{2,l) = YJ1H34] (4.12) 3 = 1 where, 'H^Vfr f E , 9 > . (4.13) The column vector 1 Hj tells us the contribution of a change in the j ' t h com-ponent of the contact coordinates to the relative spatial velocity of the two objects. Letting 1H be a matrix whose columns are 1 Hj, j = 1.. .5, we can write Equation 4.12 in matrix notation as 1(/>(2,l) = 1Hq. (4.14) This 6 x 5 matrix H relates coordinate velocities to relative spatial velocities. It can be transformed to any other convenient frame by left multiplying with an appropriate adjoint. Note that when H is written in one of the contact frames it takes on a simpler form. Its bottom row becomes zero because a velocity in the z direction is not achievable through a change in the contact coordinates (a non zero z velocity involves breaking contact). The upper 5 x 5 sub-matrix of H is in general dense. In a contact frame, however, H will have several zeros if one of the surfaces is flat and without twists in its equiparameter lines. If there is near conforming contact between the surfaces, H becomes poorly conditioned. When there is conforming contact the upper 5 x 5 sub-matrix becomes singular. Consider the example of a ball in a socket. If the surfaces have almost but not quite the same shape then very large parameter velocities are needed to account 26 s for small spatial velocities. If the surfaces conform exactly then we have a lower number of degrees of freedom and need a new parameterization for the contact (see [1] for an example of a cylinder on a flat surface going from conforming contact to point contact). Although the derivation shown in Equations 4.9-4.14 is quick for demonstra-tion, it lacks the convenience of the finer details necessary for an implementation. Appendix A contains a detailed derivation more suitable for implementation as it emphasizes ease of computation. 4.3 Const ra ined D y n a m i c s Let body 1 be free to move in contact with body 2 which is fixed. The Newton-Euler equations for a rigid body (see [36, 32]) can be written using spatial vectors as fext + fc = MJ>- [<t>}TM4> (4.15) where fext is the external force, fc is the constraint force, 4> is the spatial velocity of the body relative to the world (inertial) frame, M is the spatial inertia, and all coordinates are with respect to the body fixed frame, frame 1. Note that <j> is the time derivative of cf> in the body frame. We can combine this equation with the non-holonomic constraint 1<f>=1Hq (4.16) to get a differential algebraic equation. Differentiating the constraint once, in the body frame, we get 1<j>=1Hq+1Hq, (4.17) where *H = j c A d lcH + JlcAd lcH. (4.18) 27 Combining Equation 4.15 with the constraint Equation 4.17 (instead of Equation 4.16) gives us an ordinary differential equation instead of a differential algebraic equation. In the frictionless case the constraint force must not do any work on the system. The result of this is that fc will be orthogonal to all of the degrees of freedom. This can be expressed as HT fc — 0. If interpreted in the contact frame, fc is easily seen to be perpendicular to the surface as the first five rows of H are linearly independent and the sixth row is zero. We will write these equations in the following matrix form below where all quantities are expressed in frame 1 coordinates, so we drop the frame label for clarity. / / M o\ { f \ f b\ 0 V"5 I 0 H -a V 0 / (4.19) b=-[<f>]TM<f>-fpxt a = H q (4.20) (4.21) Performing block row elimination on the last row produces (I M 0 f b ^ 0 1 H y0 0 HTMH j \ 9 J —a Here, r = HT(fext + [4>]TM<f>-MHq). (4.22) The last row yields an ordinary differential equation in the contact coordi-nates which can be integrated directly; it's equivalent to the O D E q = {H1 MH)~lHT(fext + [<t>}lM(t> - MHq). (4.23) 28 Note that although HTMH is an invertible 5x5 matrix we use a matrix factorization instead of explicit inversion. To summarize, we choose our constraint as the description of how two sur-faces move when they are in contact. Because we constrain the motion of our objects to remain in contact, we will achieve this type of motion when we integrate while ap-plying an external force. In this approach we avoid stabilization since we evolve the system in a set of reduced coordinates which parameterize the constraint manifold exactly. Recall that we only wanted a unilateral non-interpenetration constraint. Al-though the contact kinematics equations are much less complicated than the distance computation required for a position level constraint, our formulation effectively gives us a bi-lateral constraint. We resolve this problem by monitoring the constraint force to check for when the surfaces are breaking contact. 4.4 F r i c t i o n Friction is a necessary component for realism in any physically-based animation. We present a simple way of incorporating dry Coulomb friction into our algorithm. Sections 4.4.2 through 4.4.4 discuss ideas related to static friction which we have not implemented. 4.4.1 Dynamic Friction With only a single point of contact between two contacting bodies and isotropic friction, we can make the assumption that the friction force opposes the relative velocity. When the relative translational velocity is non zero, we are in sliding or dynamic friction mode. When the magnitude of the relative translational velocity falls below a small threshold, we might consider that the contact enters sticking or non-slip or static friction mode. In this mode the system has fewer degrees of 29 freedom as the contact will move at the same speed across both surfaces. We discuss how no-slip friction could be added to our algorithm in Section 4.4.2. See also [17] for more on no-slip contact evolution. We observed results which are adequate for many graphics applications by allowing the surfaces to be in frictionless contact when the relative translational velocity falls below a threshold. When the bodies are close to a static friction configuration they will appear to be sticking but will in fact be slipping on one another and the contact will slowly drift from its current position. For translational friction due to sliding We have the equation, /t = - H / n | , « / 0 , (4.24) where ft is the tangent translational force due to friction, / „ is the normal force, and v is the normalized relative translational velocity. Summing tangent and normal components gives us the translational component of the constraint force. This equation becomes very simple when written in coordinates of a con-tact frame since the constraint wrench in coordinates of a contact frame is already decomposed into a component normal to the constraint manifold and a tangential component due to friction. The normal component is simply the z translational component (the 6th entry of i c f c ) . The remaining components are due to friction. Let lcfCt be the i t h component of the constraint wrench written in frame lc. Equation 4.24 becomes J7c4 = vx lcfc6, (4-25) lcfC5 =-uvy lcfC6. (4.26) We may also like to have friction slow down an object which is rolling or spinning. This is most commonly achieved through viscous damping due to air resistance. It is not difficult, however, to have a rolling or spinning friction torque proportional to the normal force instead of the angular velocity (there is physical justification for this in some special cases, see [31]). 30 When written in the contact frame, the relative angular velocity, u, decom-poses into a rolling component, p =f (ux1ujy)T, and a spinning component, uz. Letting p,r and Ms represent the coefficients of rolling and spinning friction we can write, C r ic r Jc\ —; — Mr Px Jce,i fc2 = — Mr Py fcei 7c3 = ~Ms sgn(wz) lcfC6. (4.27) (4.28) (4.29) Recall that the constraint wrench in Equation 4.15 is in coordinates of frame 1 . Thus we will need to left multiply by the appropriate adjoint. We can use the adjoint and Equations 4.25-4.29 to write our new constraint equation in matrix form. def Mr Px Mr Py Hs sgn(w^) M vx \ V Vy ) (4.30) '5x5 0 \MT X = - o 5x5 / i ; c A d T h (4.31) The left hand side extracts the tangential component of the constraint force, while the right hand side computes what the frictional component should be based on the normal restoring force. There are only five rows in the matrix because we do not place any restrictions on the normal force. The normal force is determined solely by the fact that the surfaces are not allowed to interpenetrate. The essence of this equation is to make the constraint wrench point in a slightly off normal direction such that it opposes the relative motion of the bodies. If the normal component 31 needs to be larger to prevent the bodies from interpenetrating the tangential portion will also become larger due to the larger normal force. Bringing everything to the left hand side in Equation 4.31 we get, F ' f c - 0, (4.32) , ,-, def where b = ^ 5 x 5 n LAdT. (4.33) V ••) The 5 x 6 matrix F thus replaces HT in Equation 4 .19 . Repeating the block row simplification performed in Section 4.3 on the modified matrix results in the following equation which we use instead of Equation 4 .23 . FMHq = F(fext + [4>]TM(p - MHq) (4.34) Note that F is a function of <p. In Equation 4 .34, a time explicit integration method would use the relative spatial velocity from the previous time step. This is what we have implemented. 4.4.2 No-Slip Friction No-slip friction imposes a non-holonomic unilateral constraint on motion. Because it is non-holonomic, it only places a restriction on the contact coordinate velocities. We still need the five dimensional contact coordinates, however, to keep track of the configuration of the contacting objects. The contact coordinate velocities on the other hand are no longer independent since they only have three degrees of freedom. In this section, unless otherwise specified, we let H be in a contact frame (say lcH, in the contact frame lc). Let i f i . . ^ . . / be the sub-matrix consisting of rows i through j and columns k through /. For example, i?4..5,i..2 is the portion of H which relates the contribution of the contact coordinate velocities (s, i) to the tangential translational velocity of the contact, (vx,vy). 32 When the two objects in contact are moving without slipping, the relative translational velocity of the bodies in the contact frame will be zero. We can use this fact to find a relationship between (s, t) and (u,v). Because ip does not contribute at all to this velocity we can write this relationship as H. 4..5,1..2 -H. 4..5,3.-4 (4.35) (4.36) where D =F - H4}5 1__2H4..5,3.A- (4.37) We can now rebuild a smaller version of H which relates the relative spatial velocity of the objects to a set of independent contact coordinate velocities. That is, we can write u \ I u = #1..6,1..2 D | + #l..6,3..4 M + ffl..6,5 V>, V I \ V , which has the matrix form = H' (u\ V (4.38) (4.39) where H' def (^1..6,1..2 D + ffi..6,3..4) Hl..6,5 I • (4-40) V i ) The no-slip equations of motion are now simply the equations in Section 4.3 with all occurances of H and q replaced with this 6 x 3 matrix H' and the reduced 33 set of contact coordinates velocities (ii, u, tp)T from Equation 4.39. We solve this smaller system to compute the reduced contact coordinate accelerations. We still need to integrate the original q, however, to get the new contact configuration q. This is done at each derivative computation step by computing q using Equation 4.36. We use the reduced contact coordinate velocities resulting from the integration of the reduced contact coordinate accelerations. 4.4.3 No-Slip and No-Spin Friction We can likewise additionally restrict the free object from spinning. This reduces the degrees of freedom to two. We get no-spin by requiring u>z to be zero in the contact frame. We can write this requirement as #3,1.-2 Q + #3,3.-4 + # 3 , 5 i> = 0. (4.41) Since lcH3<5 = 1 (see Section 4.2 or Appendix A for more details), and using D from Equation 4.36 we can write i> = - ( f f 3 , i . . 2 D + H3t3.A) . (4.42) Thus Equation 4.39 with the additional requirement of no spin becomes, (4.43) where H" 1,11 #1..6,1. .2 D + i?l . .6 ,3 . .4 + Hi..6>5(-(H3ii„2 D + H3t3.A) V : / This 6 x 2 matrix along with u and v replace H and q in our dynamics equations of Section 4.3. After integrating the computed u and v to get u and v, we use Equations 4.36 and 4.42 to compute q. 34 4.4.4 Changes in Friction We can simulate a better friction model by changing between our static and dynamic friction models at appropriate times. As mentioned before, although our sliding fric-tion model can do a reasonable job imitating no-slip friction, objects which appear to be at rest will be slowly moving. To switch between models we must monitor the constraint force and slip velocity. When in no-slip friction mode we check to see if r f i < Static (4-45) I Jn\ Recall that the tangential and normal components are easily extracted from the constraint force when it is written in the contact frame. If the inequality is not satisfied then we switch to sliding friction. When in sliding friction mode, if the magnitude of the slip velocity falls below some threshold then we switch back to no-slip friction. Note that when switching to slipping friction, the objects should be given a chance to start sliding before switching back to no-slip friction. The no-spin condition suggests that the contact area is large enough to trans-mit frictional torques. In this case we might also want to compare the torque about the z axis with the magnitude of the normal force. That is, if the inequality • f ? < A^ spin static (4.46) \Jn\ is violated then we switch to sliding friction. 4.5 Traversing Pa t ch Boundar ies Object models are commonly a collection of surface patches. In Figure 4.3, the bowl at the left consists of 56 triangular patches and the plate on the right consists of 192 triangular patches. A robust method of evolving a contact across patch boundaries is needed. 35 Figure 4.3: A bowl shape and a complex shallow bowl showing patch boundaries. Although the surfaces depicted in Figure 4.3 are subdivision surfaces, recall that any surface that can be evaluated parametrically is suitable for our contact evolution technique. For our boundary traversal method, an object can consist of any mixture of patches, from simple Bezier or N U R B S patches (see [39, 45]) to more complicated surfaces, as long as tangent plane continuity exists between all pairs of neighbouring patches. We assume that the domain of a patch is a polygon, which we describe by a counter clockwise list of points in the plane, Pi £ R 2 , i = 1.. . n. For every edge in the domain polygon we need to know both the adjacent patch and the corresponding domain edge. That is, for a given patch A and its domain edge Pf~P{\.i mapping to a curve in R 3 , we have a corresponding patch B with domain edge P?Pf+i which maps to the same curve in R 3 . We call the edge PtAP(\.i edge i of patch A. We also need to know how to find points on edge j of patch B which map to the same location in R 3 as points on edge i of patch A. One way of describing the relationship between these two boundaries is to write the boundary curve as a single parameter function. Edge i of patch A whose shape is described by cA(u, v) can have its boundary curve written as c?(u) cA(PtA+lu + PA(l-u)). (4.47) 36 Figure 4.4: Domain boundaries and functions involved in a transition Then we need to know the invertible reparameterization function g : R -» R (4.48) such that cf(u) = cf(g(u)). (4.49) In many cases this function will turn out to be quite simple. In our case, all boundary curves are of the same form and are equivalent to quartic Bezier curves. The adjacent curves are identical (the control points in the Bezier basis will be the same) giving us a very simple function g{u) = 1 - u. It is not the identity because the domain polygons are described in counter clockwise order. Thus the curves, although identical, start at opposite ends. Trimmed N U R B S surfaces, as created by many commercial design packages, would cause a problem here. Trimmed N U R B S are usually the result of a boolean 37 constructive solid geometry operation and result in sharp edges. The sharp edged curves are often described in a piecewise linear fashion in the domain since surface surface intersections result in very high degree curves making it impractical to define the boundary by a curve in the domain (see [48]). That said, we do not address curve-curve or surface-curve contacts because they require a different contact coordinate parameterization (see [1]), but [35] addresses the issue of finding boundary points in the domains of adjacent trimmed NURBS surfaces. Here are the steps for doing a domain edge traversal: • compute the state of the system at the time of crossing; • compute the contact coordinates of the system using the adjacent patch; • compute the new contact coordinate velocities. We use the following notational convention. We assume, without loss of generality, that the contact moves from patch A to patch B on body 1 while the contact remains in patch C of body 2. When the state of the system is at a patch boundary we use a superscript to denote the patch whose coordinate system is being used. For example, qA would denote the system configuration using patch A. Likewise HA denotes the contact kinematics equations using patch A. 4.5.1 Time of Crossing A boundary crossing is detected when the contact coordinates fall outside of the domain of the patch. This is easy to check when the domain polygon is convex in which case we check if the point falls to the left or right of the domain edges. If the point is on the right of any edge of the convex domain, then it is outside. Once it is determined that the contact is leaving the domain of a patch we must back up the state of the system to the time of the traversal. One approach is to revert to the previous system state and integrate with smaller time steps until the contact position in the domain of the patch being tra-38 versed is sufficiently close to the boundary. The point can then be projected onto the boundary. An alternative is to assume that the evolution of the contact is linear on the scale of the step size of the integration technique. In this case we can compute a line segment intersection to compute the new state of the system. Suppose the line segment PAP{+i is the domain boundary edge being crossed and (u,v)t(u,v)t+At is the line segment described by the inside and outside contact locations. Let a, 3 G [0,1] describe the location of the intersection on each segment as a proportion of the distance along each segment (see Figure 4.5). We can take the state of the system at the boundary to be qA = vt+pAt « A V + A * + ( i - (4-50) iA = qf+pAt « 0tf+*t + ( i - m t (4.5i) and the location of the contact on the boundary will be cf(a). Note that we define qA and qA to be the state of the system at the boundary for the convenience of not specifying the time of crossing in a subscript. As mentioned before, the superscript denotes that the contact coordinate and contact coordinate velocity use patch A. 4.5.2 Contact Location on Adjacent Patch We can write the state of the system at the boundary of patch B by using the function g defined in Equation 4.48. 'pf+1g(a) + P?(l-g(a)^ aB = 4A (4-52) V / Note that the third and fourth component of the contact coordinate remain unchanged. Note that if there is a simultaneous patch transition on the other object it will be handled after this transition is complete. Knowing the coordinates on the adjacent patch, we still need to compute tpB• This is done by adding to ipA 39 Figure 4.5: The boundary crossing as the location of intersection. a correction to reflect the angle between the x axes of the contact frames at the crossing position of the adjacent patches. Since cos~1(x'4 • xB) is computed as an angle in the range [0, TT}1, we check to see if xA X xB points in the opposite direction of the surface normal in which case we negate the correction angle. Thus we can write ^ = ^ + sgn((x^ x x s ) - z ) cos~1{xA -xB). (4.53) Note that if x'4 x x B = 0 then cos _ 1(x' 4 -xB) is either 0 or n making the sign correction irrelevant. Figure 4.6 shows the case where xA X xB points in the same direction as the normal. 4.5.3 Computing New Contact Coordinate Velocities Now that the contact coordinates q are known on both sides of the boundary we can compute qB from qA. We equate the relative spatial velocity on each side of the boundary using 2cHA and 2cHB. The spatial velocity is computed in a contact le.g. by Java 's Math.acosQ. 40 Figure 4.6: To compute contact coordinates on patch B, the angle tp has 9 added to it. We add positive 9 because xA x xB points in the same direction as the normal. frame allowing us to drop the z component of the velocity (which is zero) resulting in a system of equations which are not over-constrained. The contact frame of body 2 is used as it is unaltered while the contact frame of body 1 changes from patch A to patch B. Hence, we solve for qB in 2 c f i A . A = 2cfiB.B^ ( 4 5 4 ) where 2°H is a 5 X 5 matrix equivalent to the matrix 2cH with the bottom (zero) row removed. Knowing the system state after the boundary traversal (qB, qB), we can run the integrator again with a step size of (1 — 8) At to bring us to the target time, t + At. If additional boundary crossings occur in this time interval, then they are treated in exactly the same manner. 4.5.4 Corner Transitions We might ask what happens if the contact goes directly through a corner. In this case, although the contact may not end up in one of the adjacent patches (one sharing a boundary rather than just a vertex with our current patch), we still transfer the contact to an adjacent patch. Once in the adjacent patch the contact may evolve 41 path of contact position Figure 4.7: Transition through a corner through the corner again. The process is repeated until we end up in a patch where the contact coordinate velocity evolves the contact into the interior of the domain. If we have n patches meeting at a point, we may need to perform n — 1 traversals before we end up in the final destination patch (see Figure 4.7). This is because we take the first domain edge which intersects the state trajectory line segment. That is if the contact goes through then both PiPi+i and P; +i.P;+2 w ' l l intersect with the path. Consistently choosing the edge which ends rather than begins at the traversal point will prevent oscillating transitions between two adjacent patches. 4.5.5 Avoiding Traversal Problems To increase robustness, we also allow the contact to remain on a patch if it is outside but within epsilon of the patch boundary. This is at a small expense of accuracy, however, as the patch shape outside the boundary does not describe the shape of the surface. This epsilon extended boundary eliminates the problem of oscillating patch transitions which might occur when a contact exactly follows a patch boundary. 42 To perform transitions when the patch is crossing the epsilon extended boundary we project the contact location back onto the boundary edge before making the transition and then displace the contact location by an equivalent amount into the adjacent patch after the traversal. In doing this we must be careful not to compute a contact location in the new patch which is outside the epsilon extended boundary. 4.6 A l g o r i t h m Summary At this point it may be useful to review the entire algorithm. Everything is centered around the computation of q. For a single contact there is a constant amount of work to do this computation. While the integrator has not achieved the target time for the next frame in the animation we try to step the integrator (note that this may result in crossing a patch boundary). The integrator gives us the current state and we must compute the derivatives of this state. That is we get (q, q) and need to produce (q,q)- This only involves computing q as we can use the q provided with our current state. The whole process can be summarized briefly by the steps that follow. • First we evaluate the surfaces. That is, we must compute the position of the contact point on each surface and all the partial derivatives of order up to three. Because the basis functions of parametric surfaces are smooth and due to the symmetry of mixed partials (cjUV = c ) t m), only ten sets of basis functions are needed to do this evaluation. • Next we must compute the position and orientation of the bodies in world co-ordinates. This is necessary for computing any external forces such as gravity. We do this by computing f E from Equation 4.9 using the surface evaluations we just computed. The absolute position and orientation of the objects are also used to render the animation. 43 We then compute H with the method described in Appendix A . We transform H to the free body reference frame and compute <f> using q from the current state. We now have everything we need to solve for q. We do this by solving Equation 4.23 for the frictionless case, or by solving Equation 4.34 for the friction case. This is accomplished with an L U decomposition of HTMH in the frictionless case, or FMH in the friction case. Recall that since contact coordinates are used the constraints are automatically satisfied and there is no need for constraint stabilization. Back substituting q into the second row of Equation 4.22, we get the contact wrench fc which can be monitored to detect contact breaking. If the integrator computes a new contact coordinate which is outside the cur-rent domain boundary then we backup to the time of the crossing, traverse the boundary as described in Section 4.5 and repeat the entire process from the first step until we have arrived at the target time. 44 Chapter 5 Surface Representation We primarily use Loop subdivision surfaces to describe the boundaries of objects because of their popularity in the graphics community. Note that although subdivi-sion surfaces are similar in many ways to other popular surface representations such as NURBS, they pose some special difficulties with respect to our contact evolution technique. In this chapter we describe the considerations that are needed to use Loop subdivision surfaces with our technique. We end the chapter with a general description of the features of our implementation. 5.1 S u b d i v i s i o n S u r f a c e s Subdivision surfaces have become popular in recent years for modelling objects [12]. This is because they provide an easy mechanism for converting meshes with arbitrary topology into smooth surfaces. The original mesh provides an easy way for the user to edit the shape of the smooth surface. Another advantage is that subdivision surfaces are not hard to implement. We will give a short introduction to subdivision surfaces in this section. The SIGGRAPH course notes [12] discuss the details in much greater depth (see also Schwitzer's PhD thesis, [47]). A subdivision surface algorithm consists of rules for subdividing a polyhedral mesh where all faces have the same number of sides. The rules are applied recursively 45 to get finer and finer meshes. In the limit the rules define a smooth surface with nice continuity properties (they typically are C 1 or C2 at patch boundaries, but C°° almost everywhere). The mesh which defines the surface is usually but not necessarily a closed manifold. In the case of the Loop subdivision rules the faces are triangular, while for the Catmull Clark rules the faces are rectangular. These are the two most common types of subdivision. There are others with various interesting properties, such as the Butterfly scheme which interpolates mesh vertices. The subdivision rules define a finer mesh by moving the existing ver-tices (or not in an interpolating scheme) and inserting new vertices to make new 3 smaller faces (see figure 5.2). The rules 8 state where a vertex moves, or where a new vertex should be created, based on the positions of vertices in a small neigh-bourhood. In the refined mesh, the new vertices are called odd vertices while the Figure 5.1: Loop subdivision rules vertices which were in the original mesh are called even vertices. This naming con-vention comes from one dimensional subdivision curves because when the vertices are stored in a zero indexed array they occur at even and odd indices. Figure 5.1 shows the affine combinations, also called vertex maps, for the Loop rules. The odd vertices have their position defined by an affine combination of their position and the positions of their neighbours. The even vertices are an affine combination of the surrounding vertices. For the Loop subdivision rules, any vertex in the mesh which does not have a valence of six is referred to as an extraordinary vertex. Although there are slightly different subdivision rules at extraordinary vertices, the new position of the vertex 46 Figure 5.2: Three levels of subdivision of an octahedron. Note how the number of extraordinary vertices does not grow with subdivision. is still an affine combination of its position and the positions of its neighbours. Al-though there are multiple ways of choosing the weights (see [54]) the actual values are not important to us as we only require that the surface be smooth. Unfortu-nately, the resulting surfaces are typically less smooth at extraordinary vertices. Most meshes will have extraordinary vertices since only shapes which are topologicals a torus can have all mesh vertices with valence six. This can be seen by evaluating the Euler characteristic to compute the genus. This relationship is given in [18] by, V-E + F = 2-2g (5.1) where V" is the number of vertices, E is the number of edges, F is the number of faces and g is the genus. If a mesh consists of triangular faces and vertices of degree six, then we can count the edges and vertices from the number of faces. Each face has three vertices, each of which are shared among six faces (V = §-F). Each face also has three edges which are shared with one other face [E = §F).' Substituting these equalities into Equation 5.1 gives g equal to one, independent of the number of faces. It is interesting to note that the number of extraordinary vertices does not change as the mesh is subdivided. This can be seen in the sample subdivision of an 47 octahedron in Figure 5.2. 5.1.1 Parametric Evaluation of Subdivision Surfaces Although subdivision surfaces were not designed with parametric evaluation in mind, it turns out that they can be evaluated parametrically. Loop and Catmull Clark surfaces happen to be quite simple away from extraordinary vertices. Away from extraordinary vertices, Loop surfaces can be described as quartic box splines, while Catmull Clark surfaces can be described by bi-cubic Bezier surfaces. Because the subdivision rule for a vertex is a simple affine combination of its neighbouring vertices, a subdivision matrix can be used to describe how a neigh-bourhood of points moves at each subdivision step. In [52, 51], Stam shows that the eigenstructure of the subdivision matrix is useful for evaluating the surface next to an extraordinary vertex. The idea is to subdivide the surface until the point being evaluated is no longer in a patch with an extraordinary vertex. Subdivision to the desired depth is achieved by computing a power of the subdivision matrix. This is easily performed with the eigenvalue decomposition. In Stam's formulation each face can only have one extraordinary vertex. This is easy to guarantee as we can perform one subdivision on our mesh to get a new mesh which has this property. In the new mesh the odd vertices will have degree six and no two even vertices will be adjacent. For transitions between patches, we define the domain of a Loop patch to be the triangle with vertices Pi = (0, 0), P2 = (1, 0) and F3 = (0,1) (we treat extraor-dinary patches differently in Section 5.2). Although planar domain coordinates can be converted to barycentric coordinates for our quartic box spline representation, we do not need to do this in our implementation as we write the basis functions solely in terms of the two coordinates by simple substitution. That is, although the Loop basis functions are given by Stam (see [51]) in terms of barycentric coordinates u, v, and w, we replace occurances of w with 1 — u — v. 48 5.1.2 Non-regular Parameterization Because the contact coordinates use the parameterization of the surface, the algo-rithm is sensitive to poor parameterizations. Excessive twisting in a surface causes problems which are not always overcome with smaller step sizes. For an example where smaller step sizes work, consider a planar surface defined by a Bezier surface where the equiparameter lines twist left and right in the plane. For an object, such as a sphere, rolling in a straight line on this plane the contact coordinates will need to follow a non linear path in the domain of the twisty plane. Loop surfaces near extraordinary vertices using the natural parameterization are unpleasant. When the valence of the extraordinary vertex is less than six the partial derivatives of the patch at the extraordinary vertex are zero. That is, lim c^<6(u,v)=0, cd;s<6(0,0) = 0, (5.2) (u,v)-*(0,0) lim cdtes<6(«,«)=0, c d e g < 6(0,0) = 0. (5.3) (u,v)-*(0,0) ' When the valence is greater than six the partial derivatives at the extraordinary vertex are undefined. In this case, the magnitude of the partial derivatives become large without bound as we approach the extraordinary vertex. That is, lim |c d e s > 6(u,u)| =oo, c d e g > 6(0,0) = undefined, (5.4) (u,t;)-+(0,0) ' lim \cdtes>6{u,v)\ =oo, c d e g > 6(0,0) = undefined. (5.5) (u,i>)-»(0,0) ' ' An additional concern with both cases is that the isoparametric lines which come close to the extraordinary vertex make very sudden and sharp bends. These effects are visible in Figures 5.3 and 5.4. The patches shown in these figures are flat with the extraordinary vertices in their lower right corners. The con-trol points are placed such that n copies could be placed sharing and completely sur-rounding a degree n extraordinary vertex. The equiparameter lines become tightly bunched together near the extraordinary vertex in the degree three patch, while they become spaced further apart in the degree ten patch. The sharp bends within 49 Figure 5.3: A planar valence 3 patch showing tight bunching of equiparame-ter lines near the extraordinary vertex. Figure 5.4: Valence 10 patch showing larger spacings between equiparameter lines near the extraordinary vertex. the plane are also visible. In the valence three patch the isoparametric lines bend towards the extraordinary vertex, while they bend away from the extraordinary vertex in the degree ten case. If a contact is moving at constant speed near an extraordinary vertex, this will result in either very large or very small contact coordinate velocities (depending on the valence). The more noticeable problem, however, is that tjj will also become very large to accommodate the twists of the isoparametric lines. When a contact is evolving near an extraordinary vertex, the contact coordinate tp will have large accelerations. These accelerations may not integrate out leaving one body spinning like a top on the other. Another consequence is that H becomes poorly conditioned. The interior of a regular Loop patch, being a polynomial surface, is C°°. Note that there is normally a discontinuity in the third derivative at patch bound-aries, except when the adjacent patches describe the same function. Because of the subdivisions necessary to get a regular patch near an extraordinary vertex, there exists a spider web of discontinuities in the third derivative. One possible remedy is to reparameterize the extraordinary patch to make it regular. Stam suggests reparamaterizing based on the characteristic functions which correspond to the two partial derivatives (see the updated version of [51] which ap-pears as a part of [12]). This unfortunately involves inverting these functions which 50 is impractical since they are quartic polynomials in two variables. Additionally, derivatives of the inverted characteristic functions are needed for computing the derivatives of our reparameterized surface. Another possible solution is to use a different subdivision rule at extraordi-nary vertices. Loop recently developed a method of using a vertex map with larger support which results in a surface with bounded curvature. This research can be found in the technical report [26]. To use Loop's result we would need the eigen-structure of the new subdivision matrices. Unfortunately, the rules are identical to the original rules when the valence is three. Thus the new scheme does not alleviate the non-regular parameterization problem for the valence three case. Additionly, we did not consider this as a possible solution to our problem since this report did not come to our attention until after we implemented the solution outlined in the following sections. 5.2 Surface Replacement Given all the problems discussed in the previous section, we feel that the best option for avoiding these problems is to replace the portions of the subdivision surface near extraordinary vertices with surfaces which are better behaved. Part of our justification for replacing the surface is that subdivision surfaces are not an exact model of anything. Their main purpose is to produce a smooth surface which approximates a desired shape, using a control mesh. If we are using a subdivision surface to model a real world object, it will be at best an approximation of this object. In this line of reasoning it is not unreasonable to change the shape of the surface slightly to create a surface with better properties. This is a similar idea to what Peters does in [37]. He shows a method for building a minimal number of maximally sized N U R B S surfaces which match a Catmull Clark subdivision surface exactly, except near extraordinary vertices where they are within a given small epsilon. This would be useful if we were using Catmull 51 Figure 5.5: An four sided S-patch fit to an extra-ordinary vertex of a bowl. Note the gentle twisting in the equiparameter lines. Clark subdivision surfaces. Since we are using Loop subdivision surfaces, our solution involves cutting a hole in the surface around an extraordinary vertex and replacing it with an S-patch. S-patches, developed by Loop and DeRose in [27], are a multi sided generalization of a Bezier surface. We will give a description of S-patches in the next section. Figure 5.5 shows the result of replacing the surface near one of the degree four extraordinary vertices in the bowl on the left of Figure 4.3. Although it is hard to tell from the figure, the four sided S-patch matches the original surface quite closely being only slightly flatter than the original surface. 5.2.1 S-patches S-patches are defined by two maps. The first map is an embedding E from a polygon in the plane into an n-dimensional simplex. The second map, B, is the standard multidimensional Bezier map from the n-dimensional simplex into R 3 . When the 52 Figure 5.6: S-patch maps maps are composed, the result is a surface given by the image of the two dimensional manifold within the simplex defined by the embedding. In other words, the S-patch function mapping the polygon P to K 3 is given by the composition S = B o E as shown in Figure 5.6. The domain polygon is defined by the vertices of a regular n sided polygon centered at the origin and with the first point, Pi, lying on the positive u axis. We let the polygon be of unit area to help match the magnitude of derivatives of adjacent patches. This is desirable for contact evolution. We do not want our contact coordinate velocities to be dominated by any component as this will coincide with a less desirable condition number for H. The simplex is defined by the vertices A convenient way to think about these n dimensional vertices is to let their coordinates all be zero except for a one in the z t h position. The embedding is designed to map edges in the polygon to edges in the simplex (Loop and DeRose call this property edge-preserving). A walk around the polygon in the plane would map to a walk which goes in straight lines in each dimension and then returns to the starting point along the diagonal. 53 We will now give definitions of the maps. We will start with the multidimen-sional Bezier map. We use n dimensional barycentric coordinates, («i, • • .«„) e |o, lj (5.6) n where ^T^Mj = l, (5-7) i=i to write the Bernstein polynomials of dimension n and degree d. We also need a n dimensional multi-index for the degree d polynomials. It is defined as i= ( n , . . .in) e N n (5.8) with ii = d- (5-9) Note that N refers to non-negative integers including zero. We can now write the ith degree d Bernstein basis function as, Bernf(uu...un)=.(fjf[^. (5.10) Here we use the multinomial coefficient instead of a binomial coefficient. It is defined as, "] (5.U) Thus, given control points corresponding to the multi-indices, we can write the Bezier map as B(uu . ..«„) = ^ V? Berni(uu . ..«„). (5.12) i The embedding giving us the n dimensional barycentric coordinates we need 54 is described by the following functions. h(p) Here, p is a point in the polygon P. We use the superscripts u and u to denote the components of our two dimensional domain points. The function a; defines the signed area of the triangle Pj\P, + ip (all subscripts to the domain points are to be taken modulus n). Note that this area is zero when p falls on the edge P t P t + i . The function 7T; being the product of all but two adjacent a; functions will be zero around the entire boundary of the polygon except at the edges P 4 _iP 8 - and P{Pi+\. Each 7Tj- is normalized by dividing by the sum of all of the 7r,- functions to produce a function It is easily seen that the functions partition unity as the sum of the numerators is equal to the common denominator. The functions Z; provide the barycentric coordinates we need for the multidimensional Bezier function. Loop and DeRose write the embedding function with the letter L when the domain polygon is a regular n-gon. This embedding, which we can write as n L(p) = J2li(p)vr, (5.16) i=l has some useful properties. One such property is that k(Pi{l - u) + Pi+1u) = u. (5.17) This makes it very easy to do transitions from the adjacent Loop patches into an S-patch. The fact that domain edges map to Bezier curves (this is true for any embedding E which maps edges to edges) combined with the property stated in Equation 5.17 implies that the function g(u) introduced in Equation 4.48 is equal 55 (yu Pzu Ptu+1\ — det nV JDV T>V P ri ri+l 1 1 1 n Ki{p) M P ) + • ••• + Kn{p) (5.13) (5.14) (5.15) to 1 — u. The other important property is that we can easily maintain continuity between the S-patches and other adjacent polynomial patches (rational or non-rational). 5.2.2 Matching Continuity The blossomed version of a polynomial function is the unique affine symmetric func-tion which is equal to the original function when all the parameters are equal. For example the cubic function f(u) = u3 + u2 has the blossom f*(ui,U2,U3) = U1U2U3+ 1/3(U!U2 + U2U3 + UiU3). (5.18) The blossom is said to be affine because it is linear in each parameter when the others are held fixed. It is symmetric because the order of the parameters does not matter. The power of the blossomed form is that the representation of / in the Bezier basis is easy to compute using the blossom. The weights for the Bernstein basis polynomials are computed by evaluating the blossom with various combinations of the domain end points. For our cubic curve above, defined on the interval [0,1], we have Bezier control points /*(0, 0, 0), /*(1, 0, 0), f (1,1, 0) and /*(1,1,1). A proof of this can be found in [41]. The result is similar for multi-variable functions. In this case they are multi-affine functions rather than affine functions.since the parameters are no longer one dimensional. Ramshaw has written two tech reports, [41, 42] which discuss blossoms and their properties in much greater depth. Using the blossomed form of an adjacent triangular Bezier patch, Loop and DeRose show that continuity can be maintained between an S-patch and an adja-cent triangular patch. This technique (also demonstrated in [41] for maintaining continuity for curves and surfaces of the same type) involves placing the domains side by side and then evaluating the blossom with all combinations of the corners of 56 Figure 5.7: A triangular domain adjacent to a five sided domain our domain. Note that these points will be on the boundary or outside the domain of the blossomed function. To illustrate this further, suppose we have a triangular patch of degree d de-fined by the function / with blossom /*. To find the control points for an adjacent n sided S-patch, we build an equilateral triangle domain polygon, T, for our trian-gular patch such that it shares an edge with the S-patch control polygon. Figure 5.7 shows this where we match up the desired edge j on the triangular patch with edge i on the n sided S-patch. When written in the barycentric coordinates of the triangle T, let the points P; of the S-patch domain polygon be written as x±. We can then write the control points of the S-patch which exactly match our function / as d Vj-= /*( xi,...,xi,x2,...,x2,...xn,...,xn ). (5.19) «'l *2 in The collection of parameters corresponding to a control point and used to evaluate a blossomed function is referred to as an argument bag. Recall the ik for k = 1... n are the components of the multi-index i and should not be confused with our arbitrary choice of edge i on the S-patch. Setting all the control points for the S-patch in this way will give us an 57 Figure 5.8: An example of 5 triangular patches meeting at an extraordinary vertex. The middle refinements are shown shaded and the boundary of the S-patch is shown with a dotted line. adjacent patch which defines the exact same surface. If the S-patch is of equal or higher degree than the triangular patch then we will have C°° continuity at the join. In our case, we want to match continuity with n patches for a degree n extraordinary vertex. The n patches we use are the center refinements of the ex-traordinary patches sharing the extraordinary vertex (see Figure 5.8) . The center refinement not having any extraordinary vertices can be treated as a quartic trian-gular Bezier patch. Through our choice of where to fit the replacement surface we define the four sided domain boundary of all extraordinary patches. The domain polygon has points P i = (0,1), P2 = (0,|),P3 = (5,0) and P4 = (1,0). Note that (0,0) is always the domain point corresponding with the extraordinary vertex. Thus, each extraordinary patch at a given extraordinary vertex shares its P2P3 edge with an edge of the S-patch. 58 Figure 5.9: A non uniquely determined control point. This control point is at a distance of one from a boundary of both triangular patches A and B. Ck continuity is maintained by adjusting the positions of control vertices which are within a distance of k of the boundary. That is these control points correspond to argument bags with a maximum of k points which are not endpoints of the shared edge. By choosing to maintain C2 continuity between the degree four S-patch and the surrounding Loop patches, we completely determine the shape of the S-patch. Although all control points are defined by a boundary, they are not defined uniquely. Consider two adjacent edges, Pi-iPi and P;P;+i as shown in Figure 5.9. The two edges border on different patches, but both patches will want to set the control point corresponding to the argument bag Pt-, F,-, because the control point corresponding to this bag is at a distance of one from both boundaries. Our solution is to set these control points which have multiple possible set-tings on a first come first serve basis. The desired shape is preserved since the patches adjacent to our two consecutive edges are C 1, regardless which patch we use to set the control point. We can demonstrate why this happens using a property of the Bezier curve control polygon. Consider a corner vertex where two triangular patches meet an 59 5 4 Figure 5.10: The dark rectangle shows the tangent plane at vertex 2. Corner control points are shown as black circles while boundary control points are shown as empty circles. For clarity, this diagram shows the case of cubic triangular patches. The interior control points of the triangular patches are omitted for clarity. The segment of all the boundary Bezier curves involving vertex 2 lie in the tangent plane. This is also true for the curve defined by V 5 5 5 , V552, V522, V222 even though it may not actually lie on the five sided S-patch. 60 S-patch as shown in Figure 5.10. The corner vertex and the control points at a distance of one on the triangular patches will all lie in the same plane if the surface is G1 (there will be additional restrictions on their position if it is C 1). This is because the derivative of a Bezier curve at its endpoint is a multiple of the second last control point subtracted from the last control point. The multiple is equal to the order of the curve. All the control points of the S-patch at a distance of one from this corner vertex must also lie in the plane because of how we computed them. Consider the curve whose control polygon connects connecting the simplex vertex corresponding with the corner vertex to any other vertex in the simplex (vertex 2 to vertex 5 in the figure). This curve will be continuous at the join even though the curve as a whole may not lie in the S-patch. The control points will be the same as those we would compute if we were matching continuity for this curve alone. This is true for each curve corresponding to each direction u,- in the simplex. The result is that the embedding is irrelevant. A directional derivative taken at our join vertex (vertex 2 in the figure) with a direction s in the simplex will lie in the tangent plane. By definition, Evaluated at the corner vertex, each 4^- lies in the plane and thus their linear com-' aui r bination must lie in the plane too. Note that the function may not be regular if there is some combination of the partials equal to zero. This is only possible, however, with meshes having an irregular collection of obtuse angle triangles meeting at an extraordinary vertex. These meshes could be made suitable for surface replacement either through retriangulation or through moving control vertices. The only side effect of our first come first serve rule is some minor twisting of equiparameter lines within the surface. This is not unexpected and can be seen in Figure 5.5. Although it may be possible to adjust the size and shape of the adjacent triangular domains T to achieve equal blossom evaluations for these shared control vertices, we have not explored this avenue. Averaging the control points may also (5.20) 61 be a reasonable thing to do. The only exception to the first come first serve rule is when a control point corresponds to an edge. In this case we let the adjacent patch define the position. If we did not do this then the control point corresponding to bag (P;, P;, Pt, P;+i) on edge i, if set by the patch adjacent to edge P;_iPt in order to maintain C 1 continuity, would break C° continuity along the edge PtPj+i leaving us a surface with a hole in it. 5.3 M o r e on C o m p u t i n g Derivat ives and Blossoms Evaluating derivatives or blossoms of regular Loop patches is relatively easy. The terms of the basis polynomials have factors of the form uz and v3. For any given partial derivative we can simply compute the derivative or blossom of the terms and then evaluate the basis functions. For example, any time u4 appears as a factor within a term of a basis function we write u4 instead. We do this for all powers of u and v. If we are evaluating a partial derivative with respect to u then we set u4 to 4M 3, U3 to 3M2, and so on down to uO which we set to zero. When coding the basis polynomials it is important to include uO as a factor of terms which do not explicitly contain a power of u. For example, the term v3 should be written as v3*u0 if we are to have correct derivative computations. This method also works for extraordinary patches since they use the same basis functions. We should mention, however, that we had a problem with derivative computations of extraordinary Loop patches due to a small error in [51]. Stam points out that the magnitude of derivatives must be corrected to take into account a remapping of the domain. This involves a sign change if the point evaluation is performed within a triangular patch which points in the opposite direction of the original patch. This sign change is not necessary for even derivatives, however, as applications of the chain rule will provide and even number of negations. The code 62 in the paper negates the result if the order of differentiation is odd but only if the number of subdivisions necessary to put the evaluation point in a regular patch is also odd. The evaluation of blossoms is achieved in a very similar way. Each factor is replaced with a symmetric affine polynomial whose diagonalization is equal to the factor. Blossoms of quartic Loop patches have four parameters. Instead of u we now have u\, u2, u3, and u4. For the factors corresponding to powers of u we set our variables as follows. 110 = 1 U l = + U2 + U3 + U4) 4 U2 = ~(UiU2 + U2U3 + U3U4 + UiU3 + U2U4 + U\U4) U3 = ~{UiU2U3 + U\U2U4 + U\U3U4 + U2U3U4) u4 = U\U2u3u4 The factors corresponding to powers of v are set in the same manner. Computing the blossom of an extraordinary patch, however, does not give expected results. Recall that we match continuity with the center refinement of an extraordinary patch. We tried to join an edge of an S-patch with the line segment (0, |)(|, 0) within an extraordinary patch. This line segment corresponds to a quar-tic Bezier curve. The degree would be much higher for an arbitrary line segment, but since it is an edge of the center refinement we are guaranteed that it is only quartic and thus it will be possible for us to construct an adjacent quartic S-patch that matches continuity. We expected that we could match continuity with the center refinement in this manner without performing a mesh refinement to find the control points of the center refinement. This naive approach, however, does not yield the desired control points. We believe the control points are displaced due to the conversion of the extraordinary Loop patch's control points into the eigenbasis. Although the correct 63 control points can likely be computed by applying the inverse transformation, we chose to use the center refinement because it is easier to compute. 5.3.1 Derivatives of S-patches Because S-patches are defined by a rather complex multivariate rational polynomial, special care is required to compute partial derivatives correctly. We used basic automatic differentiation techniques to write this code by hand. We tested the resulting code by comparing evaluations of all the partial derivatives at random points with evaluations of the derivatives computed symbolically with Maple. For each assignment statement involved in a computation, an automatic dif-ferentiation programs creates a new variable to represent the derivative of the vari-able on the left hand side of the assignment and creates a new assignment statement that computes the derivative value for the new variable. For example, if we want a partial derivative with respect to u, the assignment t= l will cause the automatic differentiation program to create the new variable t u and the preceding statement tu=0. The new statement must precede the original because a variable might be used on both sides of an equation. Consider the statement t= t*x [ i ] . This statement results in the creation of the preceding statement t u = t u * x [ i ] + t * x u [ i ] . For our implementation we need both ti and v partial derivatives in all combinations up to order three. The computation of the 7r,(p) functions defined in Equation 5.14 is a repre-sentative example of applying the automatic differentiation technique. The block of code in Figure 5.11 shows how we compute the 7T;(p) functions. 5.4 Implementat ion Desc r ip t ion We have implemented this algorithm in Java using Java3D for some of the data structures and for displaying the animation. We have a custom scene description file format for describing a system. A 64 for ( in t i piuuu [ i ] p iuuv[ i ] = 0 . 0 p iuvv [ i ] = 0 . 0 p i v v v [ i ] = 0.0 p iuu[ i ] = 0 . 0 p i u v [ i ] = 0.0 p i v v [ i ] = 0 . 0 p i u [ i ] = 0.0; p i v C i ] = 0 . 0 ; p i [ i ] = 1.0; for ( int j = 0; i f ( j == i 0; i < N; i++ ) { 0.0; < N; j++ ) { j == (i+N-l)'/.N ) continue; / / note that these are s i m p l i f i e d due to the constant zero / / second and t h i r d order p a r t i a l der iva t ives of alpha piuuu[ i ] = piuuu[i] * a lpha[j] + 3 * p iuu[ i ] * a lphau[ j ] ; p iuuv[ i ] = piuuv[ i ] * a lpha[j] + p iuu[ i ] * alphav[j] + 2 * p i u v [ i ] * a lphau[ j ] ; p iuvv[ i ] = p iuvv[ i ] * a lpha[j ] + 2 * p iuv [ i ] * a lphav[j] + p i v v [ i ] * a lphau[ j ] ; p ivvv [ i ] = p ivvv [ i ] * a lpha[j ] + 3 * p i v v [ i ] * a l p h a v [ j ] ; p iuu[ i ] = p iuu[ i ] * a lpha[j] + p i u v [ i ] = p iuv [ i ] * a lpha[j] p i v [ i ] * a lphau[ j ] ; p i v v [ i ] = p i v v [ i ] * a lpha[j] + 2 * p i u [ i ] = p i u [ i ] * alphaCj] + p i [ i ] p i v [ i ] = p i v [ i ] * a lpha[j] + p i [ i ] p i [ i ] = p i [ i ] * alpha [j] ; 2 * p i u [ i ] * a lphauCj] ; p i u [ i ] * alphav[j] + p i v [ i ] * a lphavEj ] ; * a lphau[ j ] ; * a lphav[ j ] ; Figure 5.11: Code produced with automatic differentiation techniques for computing the 7T;(p) functions 65 scene description file contains information about the objects in a system, how they should be rendered, and their initial conditions. An example file can be found in Appendix B . We used javacc to create a parser for this file format which makes it easily altered or extended. Rigid bodies in this file can be defined by either Loop subdivision surfaces or they can have a limited portion of their boundary defined by a bi-cubic Bezier patch. The Loop subdivision surfaces are defined by meshes stored in obj format. An example mesh file can also be found in Appendix B . These obj files can be written by hand for simple objects. For larger objects, however, we find it more convenient to use a modelling package such as 3D Studio Max. Our interface is minimal as we only use it to test the underlying algorithm. Commands can be typed in a text box to change how objects are displayed, and to modify properties in the simulation such as friction. Interaction with the simulation is possible by activating a spring force which connects the center of mass of the free object to a point which the user controls with the mouse. We find this minimal form of interaction useful for testing. We also have an option for saving a selected time range of an animation as a Renderman rib file for rendering. The rendered images can be much more com-pelling because they show shadows that help convey where the contact is occurring. Unfortunately there is currently no support for shadows in Java3D. The ability to render images from the simulation to a file is also useful for creating animations. The ability to easily save Java3D images to a file has only recently become possible. 66 Chapter 6 S i m u l a t i o n R e s u l t s a n d E v a l u a t i o n We present simulation results with a focus on how much time is spent in various parts of the algorithm. This is important as we wish to show that this algorithm can be used at interactive rates. We will also describe several simulation sequences which should help demonstrate the success of our algorithm in simulating single contact evolution. 6.1 Timings Running the HotSpot Java virtual machine on a 350MHz Pentium II machine the entire computation of q takes about 1.2 ms. With this computation time the simu-lation can run at 15 frames per second without using all available C P U cycles. The bulk of the time is spent in Java3d code to draw the system. Table 6.1 shows the time necessary to evaluate various types of surface func-tions along with all the necessary partial derivatives for the algorithm. We per-formed all our tests on a 350 M H z Pentium II running Java 1.3. A l l times reported in this section are measured as an average of ten thousand computations to give the HotSpot virtual machine sufficient opportunity to optimize our code. Although 67 Table 6.1: Timings for calculating the ten partial derivatives necessary for the al-gorithm. Timings are with Hotspot optimizations, run on a 350MHz Pentium II. Times are an average of ten thousand successive computations and hence may ap-pear smaller than they should due to caching effects. surface type evaluation time Loop regular 0.2413 ms Loop extraordinary (degree 3) 0.1683 ms Loop extraordinary (degree 5) 0.1933 ms Loop extraordinary (degree 7) 0.2223 ms Loop extraordinary (degree 10) 0.2654 ms S-patch (3 sides) 0.1632 ms S-patch (4 sides) 0.4466 ms S-patch (5 sides) 1.1006 ms S-patch (10 sides) 20.0248 ms Bi-cubic Bezier 0.0801 ms timings in the table may be smaller than in practice due to caching effects, all other times reported in this section do not have this problem as they were measured during an actual simulation. We found that evaluating the surface functions takes up a substantial portion of the time spent by the simulator. In our implementation these functions have quite a bit of room for optimization. Surface representations which are quick to evaluate are thus preferred for our algorithm. For two contacting cubic Bezier surfaces, the program computes H in 0.3485 ms on our test machine. Note that Table 6.1 reveals that half of this time is spent on surface evaluations of the two patches. The routine which computes q takes about 0.235 ms (this includes the L U factorization and back substitution). The total computation time for the derivatives comes to about 0.9 ms in the Bezier on Bezier case. The time which is unaccounted for comes from checking boundaries, computing object positions, velocities, and external forces. 68 6.2 B o w l on F la t Frict ionless Surface The frictionless setting is useful for evaluating the correctness of a large part of our algorithm. Without friction, the center of mass of an object on a flat surface will not be able to have any horizontal translational motion. Observing purely vertical motion of the center of mass validates much of our implementation. We simulated the system shown in Figure 6.1 consisting of an oddly shaped bowl sliding on a flat frictionless surface. The bowl is a Loop subdivision surface while the flat surface is a Bezier patch. The bowl's center of mass is close to its base. We set gravity as the only external force in the simulation. The upper strip shows a wire frame outline of each triangular patch so that the center of mass (the large sphere) and the contact point (the small sphere) can be seen. The triangular patch which is currently in contact is drawn with greater detail. We observed that the center of mass has a purely vertical motion. This result is shown in Figure 6.2. Over time the horizontal position of the center of mass will change, however, as our numerical method has limited precision (this is not visible in the figure because of the error introduced by small cusps at S-patch boundaries). To prevent drift, we would either need to identify this special situation and formulate it as a three degree of freedom system, or we would need to use stabilization to keep the center of mass above the desired position in the plane. The simulation will not handle the case of the bowl lying upside down on the table as the level lip of the bowl would establish contact with the flat plane all at once. In the simulation above, however, the system lacks the energy to reach this configuration because of the initial conditions. 69 Figure 6.1: Simulation of a bowl sliding on a flat surface. Our system is validated by the vertical motion of the center of mass which can be seen as the larger sphere in the upper sequence. The smaller sphere shows the location of the contact point. 7 0 _ Q -| I i i i i i I 0 2 4 6 8 10 12 time Figure 6.2: Graph showing the position of the center of mass in world coordinates over time. The solid line is the z position while the dotted lines are the x and y positions. The x and y positions do not move until a transition into an S-patch. This suggests a small problem with S-patch transitions. We attribute this problem to a cusp at the boundary due to an error in our implementation. 71 Figure 6.3: A lumpy marble rolling on a plate. Frames are one second apart. The free object can be seen to slide almost to the center of the plate before it starts to roll. 6.3 B lobs in Bowls Because our simulator only handles single contact, we can only examine the small time segments involving a single continuous contact for interactions between arbi-trary objects. An alternative is to restrict ourselves to objects which remain in contact and can only ever contact at a single point. This is an interesting test case because we can run our simulator for an arbitrary length of time. Although these types of system may have only limited interest, they are more interesting than the curved-flat surface interaction in Section 6.2. They help validate our implementation as the matrix H is more complex in this case. The matrix in the flat surface case has more zero entries because the contribution of a flat surface's parameter velocities to u is zero, except for LOZ which may be non-zero if the equiparameter lines of the flat surface curve within the plane. Figure 6.3 shows an example of curved-curved contact between a lumpy mar-ble and a shallow bowl. Both surfaces are Loop subdivision surfaces. Again gravity is the only external force, however due to friction the marble quickly starts to roll 72 Figure 6.4: Two different rattleback tops with asymmetric shape. after a short period of sliding towards the center of the bowl. 6 . 4 Ra t t leback S imula t ion Spinning tops, in most cases, maintain a single continuous contact with the surface on which they are spinning. This makes tops an ideal object for us to simulate with our algorithm. Rattleback tops, also known as celts or wobblestones, are interesting because they can reverse their spin. Some rattlebacks will reverse their spin in both directions while others have a spin bias and will only reverse their direction if spun in the direction opposite to their preferred direction. As shown in Figure 6.4, rattlebacks have a long elliptical shape. Those in the figure have an asymmetric shape, but this is not necessary for spin reversal to occur. When a rattleback is spun in the appropriate direction, its spin will slow and it will begin to wobble. The wobbling becomes larger until the top is no longer spinning and is instead rocking back and forth. The top then starts to spin in the 73 tps-V; in ' 'W J " other direction and spins faster as the wobbles become smaller. It is due to the shape and/or inertial asymmetry of the top that the spinning energy is converted to a rocking motion and then back to a spinning motion. Rattlebacks have been studied for a long time. The original paper which analyzed the motion of a rattleback appeared in 1896. More recently, Garcia and Hubbard in [16] explain how spin reversal can occur in one or both directions and discuss the limitations of the models in previous investigations. The limitations of previous work are mostly due to assumptions such as no-slip friction, assumed shape, and poor modelling of energy dissipation. For our simulation we build a very sim-ple rattleback model. Our model shown in Figure 6.4 consists of a single bi-cubic Bezier patch. It is longest in a direction 10 degrees off of the x axis and its width is one third of the length. Its curved shape comes from ele-vating the outer control points. We spin the model about its z axis. A scene description file ^ . „ _ ~ , , , , , Figure D . 5 : Our rattleback model, including our rattleback model can be found top and front views in Appendix B. Although the top is symmetric, we set the inertia tensor such that it is not in alignment with the planes of symmetry. In a physical model this can be achieved through a nonuniform mass density. The reversal effect is due to the misalignment of two of the principle axes of inertia with the principal axes of curvature at the point of contact. With this simple model, using sliding friction and rotational dissipation we can simulate single or multiple spin reversals depending on the coefficients of rota-tional dissipation. The preferred direction of spin can be found by giving it a tap on one end to 74 Figure 6.6: Tapping a rattleback to start it spinning. Previous positions are shown in grey. The top quickly starts spinning in a clockwise direction. Figure 6.7: Spin start of a rattleback. Previous positions are shown in grey. The direction of the spin is initially counter clockwise and quickly changes to clockwise. start it rocking. Even rattlebacks which reverse in both directions have a preferred spin direction. They take longer to reverse their spin when spinning in the preferred direction. Figure 6.6 shows a sequence of frames from a tap start, while Figure 6.7 shows a sequence of frames where it is spun to start. Several previous positions of the top are drawn in grey to give an indication of the direction of motion. In the latter simulation we start the spin slightly off the center of the top. This is necessary because a perfect spin about the z axis results in a stable rotation. These figures have coefficients of friction set as follows. pi = 0.3 Mr = 0.0 ns = 0.0 If we set pir = 0.003 and ps = 0.0001 we can prevent the second spin reversal. If all friction coefficients are set to zero we do not get any reversal at all. 75 Figure 6.8: A sequence showing a second spin reversal. Previous positions are shown in grey. The top spins clockwise at 23 seconds and counter clockwise at 40 seconds. 76 Chapter 7 Conclusions We have described a fast contact evolution algorithm for single contact between piecewise parametric surfaces. We first derived contact kinematics equations for arbitrary regular parametric surfaces. Then we showed that using these equations as acceleration constraints we can formulate the contact dynamics as an implicit O D E in the contact coordinates. The resulting O D E can be easily integrated using explicit integrators, without the need for constraint stabilization. We can easily incorporate a Coulomb friction model into our formulation. Although our implementation of friction only provides sliding friction we find that we can still achieve a good approximation of no slip friction effects such as pure rolling. We also show how our formulation can be modified for no-slip and no-spin friction. Because objects are made up of multiple patches, we derive a method for evolving a contact across patch boundaries. We also address special transitions such as going through a corner and propose extending patch boundaries by a small amount to increase stability through prevention of oscillating transitions. We implement our algorithm using Loop subdivision surfaces. Unfortunately, the natural parameterization of Loop subdivision surfaces near extraordinary ver-tices is not regular. Our solution to this problem is to replace the portions of the 77 subdivision surface near extraordinary vertices with an n sided patch called an S-patch. This involves computing control points of the S-patch so that continuity is maintained at the join. Several example simulations demonstrate the success of our formulation. We use a frictionless system to partially validate our implementation, and we use a simulation of a marble in a bowl to test our friction implementation. The results also show that our technique predicts the reversal behaviour of a rattleback top. We also measured the time used performing various calculations for our sim-ulation and discovered that a large proportion of the time, sometimes more than half, is spent computing partial derivatives of the parametric surfaces. This sug-gests that surfaces which are less expensive to evaluate are more desirable with our technique. The main limitation of our algorithm is that we only treat single contact. Although only a small piece of the puzzle, we feel this method of dealing with smooth surface contact evolution is both an interesting and useful contribution to the rigid body dynamics community. 7.1 Future W o r k Extending our implementation to handle other types of contact would make it much more useful. Flat conforming contacts as well as curve-surface and curve-curve contacts occur frequently in reality. Likewise, multiple contacts between two rigid bodies is a difficult problem which also occurs frequently and should be investigated. Lastly, an implementation capable of simulating many different kinds of contact would greatly benefit from a general framework for changing contact types. 78 B i b l i o g r a p h y [1] M . Anitescu, J . Cremer, and F . Potra. Formulating 3d contact dynamics prob-lems. Mechanics of Structures and Machines, 24(4):405-437, 1996. [2] U . Ascher. Stabilization of invariants of discretized differential systems. Nu-merical Algorithms, 14:1-23, 1997. [3] U . Ascher, H . Chin, L . Petzold, and S. Reich. Stabilization of constrained me-chanical systems with daes and invariant manifolds, the Journal of Mechanics of Structures and Machines, 23(2):135-157, 1995. [4] U . Ascher, D. K . Pai, and B . Cloutier. Forward dynamics, elimination methods, and formulation stiffness in robot simulation. International Journal of Robotics Research, 16(6):749-758, December 1997. [5] D . Baraff. Analytical methods for dynamic simulation of non-penetrating rigid bodies. In Proceedings of SIGGRAPH '89 (Boston), volume 23 of Computer Graphics Proceedings, Annual Conference Series, pages 223-232. A C M SIG-G R A P H , July 1989. [6] D . Baraff. Curved surfaces and coherence for nOn-penetrating rigid body sim-ulation. Computer Graphics, 24(4):19-28, August 1990. [7] D . Baraff. Issues in computing contact forces for non-penetrating rigid bodies. Algorithmica, 10:292-352, 1993. [8] D . Baraff. Fast contact force computation for nonpenetrating rigid bodies. In A . Glassner, editor, Proceedings of SIGGRAPH '94 (Orlando, Florida, July 24~ 29, 1994), Computer Graphics Proceedings, Annual Conference Series, pages 23-34. A C M S I G G R A P H , A C M Press, July 1994. ISBN 0-89791-667-0. [9] J . Baumgarte. Stabilization of constraints and integrals of motion in dynamical systems. Computer Methods in Applied Mechanics and Engineering, 1:1-16, 1972. 79 [10] C. Cai and B. Roth. On the spatial motion of rigid bodies with point contact. In IEEE Conference on Robotics and Automation, pages 686-695, 1987. [11] J . F . Cremer. An Architecture for General Purpose Physical System Simulation - Integrating Geometry, Dynamics, and Control. PhD thesis, Cornell Univesity, 1989. [12] T. DeRose. Subdivision surface course notes. S I G G R A P H Course Notes, 1998. [13] D . DiFilippo. The AHI , an integrated audio haptics interface. Master's thesis, University of British Columbia, 2000. [14] D . DiFilippo and D . K . Pai. The AHI : An audio and haptic interface for contact interactions. In ACM Symposium on User Interface Software and Technology, 2000. [15] R. Featherstone. Robot dynamics algorithms. Kluwer, 1987. [16] A . Garcia and M . Hubbard. Spin reversal of the rattleback: theory and ex-periment. Proceedings of the Royal Society. London. Series A. Mathematical, Physical and Engineering Sciences, 418(no. 1854):165-197, 1988. [17] Si Goyal. Second order kinematic constraint between two bodies rolling, twist-ing and slipping against each other while maintaining point contact. Technical Report TR89-1043, Cornell University, Computer Science Department, October 1989. [18] M . J . Greenberg and J . R. Harper. Algebraic Topology, A First Course, chapter Betti Numbers and Euler Characteristic. Number 58 in Mathematics lecture note series. Addison-Wesley, 1981. [19] J . K . Hahn. Realistic animation of rigid bodies. ACM Computer Graphics, 22(4):299-308, 1988. [20] G . Hegron. Rolling on a smooth bi-parametric surface. In Eurographics Work-shop on Animation and Simulation, pages 205-213, 1991. [21] A . Jain. Unified formulation of dynamics for serial rigid multibody systems. Journal of Guidance, Control, and Dynamics, 14(3):531-542, 1991. [22] Y . Jia and M . A . Erdmann. Observing pose and motion through contact. In Proceedings of the IEEE International Conference on Robotics and Automation, 1998. 80 [23] D . E . Johnson and E . Cohen. A framework for efficient minimum distance computations. In Proc. International Conference on Robotics and Automation, pages 3678-3684. I E E E , May 1998. [24] C. W . Kilmister and J . R. Reeve. Rational Mechanics. American Elsevier Publishing Company, Inc., 1966. Library of congress catalog card number 66-11823. [25] Z. L i and J . F . Canny. Motion of two rigid bodies with rolling constraint. IEEE Transactions on Robotics and Automation, 6(l):62-72, 1990. [26] C. T. Loop. Convex triangular subdivision surfaces with bounded curvature. Technical report, Microsoft Research, 2000. [27] C. T. Loop and T. D . DeRose. A multisided generalization of bezier surfaces. ACM Transactions on Graphics, 8(3):204-234, July 1989. [28] C. Lubich, U . Nowak, U . Pohle, and Ch. Engstler. Mexx - numerical software for the integration of constrained mechanical multibody systems. Tech. report SC92-12, Konrad-Zuse-Zentrum fur Informationstechnik, Berlin., 1992. [29] B . V . Mirtich and J . F . Canny. Impulse-based dynamic simulation of rigid bodies. In Symposium on Interactive 3D Graphics, 1995. [30] D . J . Montana. The kinematics of contact and grasp. International Journal of Robotics Research, 7(3):17-32, 1988. [31] D. F . Moore. The Friction and Lubrication of Elastomers, chapter Chapter 4, Rolling and Sliding. Pergamon Press Ltd. , 1972. [32] R. Murray, Z. L i , and S. S. Sastry. A mathematical introduction to robotic manipulation. C R C Press, 1994. [33] Ju. I. Neimark and N . A . Fufaev. Dynamics of Nonholonomic Systems. Amer-ican Mathematical Society, 1972. [34] D. D . Nelson and E . Cohen. User interaction with cad models with nonholo-nomic parametric surface constraints. International Mechanical Engineering Congress and Exposition, November 1998. Haptic Symposium. [35] D. D . Nelson, D. E . Johnson, and E . Cohen. Haptic rendering of surface-to-surface sculpted model interaction. In Proceedings of the ASME Dynamic Systems and Control Division, volume DSC-Vol . 67, pages 101-108, 1999. 81 [36] D . K . Pai, U . M . Ascher, and P. G . Kry. Forward dynamics algorithms for multibody chains and contact. In IEEE International Conference on Robotics and Automation, pages 857-863, 2000. [37] J . Peters. Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values. In SIGGRAPH, 1998. [38] F . Pfeiffer and C. Glocker. Multibody Dynamics with Unilateral Contacts. Wiley series in nonlinear science. John Wiley and Sons. Inc., 1996. [39] L. Piegl and W . Tiller. The NURBS Book. Monographs in Visual Communi-cation. Springer, second edition, 1997. [40] P. J . Rabier and W . C. Rheinboldt. On the numerical solution of the euler-lagrange equations. Technical report, University of Pittsburgh, February 1993. [41] L. Ramshaw. Blossoming: A connect-the-dots approach to splines. Technical Report 19, Digital Systems Research Center, 130 Lytton Avenue, Palo Alto, California 94301, June 1987. [42] L. Ramshaw. Blossoms are polar forms. Technical Report 34, Digital Systems Research Center, 130 Lytton Avenue, Palo Alto, California 94301, January 1989. [43] R. E . Roberson and R. Schwertassek. Dynamics of Multibody Systems. Springer-Verlag, 1988. [44] G . Rodriguez, A . Jain, and K . Kreutz-Delgado. A spatial operator algebra for manipulator modeling and control. The International Journal of Robotics Research, 10(4):371-381, 1991. [45] D. F . Rogers. An introduction to NURBS: with historical perspective. Morgan Kaufmann, 2001. [46] J . Sauer and E . Schomer. A constraint-based approiach to rigid body dynamics for virtual reality applications. In ACM Symposium on Virtual Reality Software and Technology, pages 153-161, 1998. [47] J . E . Schweitzer. Analysis and Application of Subdivision Surfaces. PhD thesis, University of Washington, Seattle, Washington 98195, August 1996. Technical Report UW-CSE-96-08-02. 82 [48] T. W. Sederberg, D . C. Anderson, and R. N . Goldman. Implicit representation of parametric curves and surfaces. Computer Vision, Graphics, and Image Processing, 28:72-84, 1984. [49] J . M . Snyder. An interactive tool for placing curved surfaces without interpen-etration. In R. Cook, editor, SIGGRAPH 95 Conference Proceedings, Annual Conference Series, pages 209-218. A C M S I G G R A P H , Addison Wesley, August 1995. held in Los Angeles, California, 06-11 August 1995. [50] J . M . Snyder, A . R. Woodbury, K . Fleischer, B . Currin, and A . H . Barr. Interval method for multi-point collisions between time-dependent curved surfaces. In SIGGRAPH 93 Conference Proceedings, pages 321-334. A C M S I G G R A P H , 1993. [51] J . Stam. Evaluation of loop subdivision surfaces. In SIGGRAPH, 1998. In-cluded on course notes C D - R O M . [52] J . Stam. Exact evaluation of catmull-clark subdivision surfaces at arbitrary parameter values. In SIGGRAPH, 1998. [53] D . Stewart and J . C. Trinkle. An implicit time-stepping scheme for rigid body dynamics with coulomb friction. In IEEE International Conference on Robotics and Automation, pages 162-169, 2000. [54] J . Warren. Subdivision methods for geometric design. Rice University, Novem-ber 1995. [55] A . Witkin and D. Baraff. Physically based modeling: Principles and practice. S I G G R A P H Course Notes, 1997. 83 Appendix A Contact Kinematics Derivation We use the same setup as described by Equations 4 . 5 - 4 . 8 in Section 4 . 2 . As the contact point moves, so do the contact frames. The spatial velocity of the contact frame lc relative to the body frame 1 in terms of contact parameter velocities, has a particularly simple form in the contact frame. If lc(f>(l,lc) = (uT,vT)T, then the spatial velocity of the contact frame lc relative to the body frame 1 in coordinates of lc is M V \ i c p 1 F lCf- 1 i E icE,s s + , E J C E , f t 0 0/ (A. l ) Let 0 be the leading 3 x 3 sub-matrix of j c E , and we'll do the matrix multiplications in Equation A . l to get a closer look at the contributions of s and i to u and v. Recall that the rightmost column of \CE is the location M = o T 0 , , s + eT©,t i, v = 0 T c , s s + eTct i ( A . 2 ) Lets first look at u> by examining the skew symmetric matrix 0 T 0 , S ( 0 T 0 , t will be very similar). / X • X . X-V X - 7 . \ 0 T 0 S = y -x , s z • x x - y , s y -y , s z - y,5 y z,s z • z (A.3) 8 4 The least complicated expressions relating OJ to s are those which are boxed in Equat ion A . 3 (x i S is less complicated than z ) S which is less complicated than y s ) . Thus we can write u as, UJ = -z • X s + /-,.It\ Now looking at the 0 T c s and © T c * components we can write, / x - c \ Q1 C,s = y-c, (x.cA 0 V 0 J and QTc t = Equat ions A . 4 and A . 5 combine to give °H\. ( lc T T def til = yz, s - y z - z x 5 - z - x , t y • x , 5 y • * t x • c s x • ctt o y c i \ o o j y • c,t V Z C J .A /x • C i \ y -c,i v ° / (A.4) (A.5) (A.6) lc4>(l,lc) = lcH1 (A.7) We can analogously define °H2. Th i s transforms to frame lc as lcH2 = iccM 2cH2 RTJ, 0 i Z c H2. (A.S) 0 R^/ Final ly , relative spatial velocity of the two contact frames is a pure rotation 85 c<f>(lc, 2c) (A.9) about the surface normal, i.e., / o \ 0 -1 0 0 We can now compute the contact matrices Hk as follows. The relative spatial velocity of the two bodies is given by 4>(1,2) = 4>{l,lc)+4>{lc,2c) + <f>{2ci2) = <f>(l, lc)+ <i>(lc,2c) - 4>{2,2c). Substituting Equations A.7, A.8 and A.9 we have / : : : \ lC jj lC JJ IC TT V ; ; : j . def lc rT . q = Hq (A.10) The contact matrix H can now be transformed to any convenient frame, for instance, frame 1: JH = 'lcAd lcH. In Section 4.3 we need H as we will take the time derivative of Equation A.10. Once lcH\ and 2°H2 are computed with the current state (q and q) using the chain and product rules, we can write ( : lcHx -(^A'd 2 C F 2 + ; cAd 2 ci7 2) 0 (A.ll) 7 86 Appendix B Example Files This appendix shows example files used with our implementation. The scene de-scription file format is the only new format that we use. Note that the subdivision surfaces use the standard obj format while stripified versions of these objects are stored in objf format as created by the program stripe. Texture files for Java3d and Renderman are in standard jpeg and tiff. Note that Renderman compliant rendering programs will generate the txt format texture files they need. B . l Scene Descr ip t ion F i le Here is an example scene description file. Note there is a description of a Loop subdivision surface commented out. // This f i l e simulates a rattleback view VI { Translate [0 -0.355 1.2]; Rotate -130 [ 1 0 0]; ribFileName = " r b . r i b " ; frameBaseName = "z:/pgkry/frames/rb"; Format = "300 225 -1"; // size of the image... what i s -1? PixelSamples = "2 2"; >; World RBSim { LightSource LI ambientlight { 8 7 i n t e n s i t y = 0.3; colour = [ 0.1 0.1 0.1 ]; index = 0; >; LightSource L2 arealight { point = [ 0 40 40 ]; po i n t i n t e n s i t y = 1.0; colour = [ 1 1 1 ] ; in t e n s i t y = 3585; polygonp = "[-10.0 10.0 40 10.0 10.0 40 10.0 -10.0 40 -10.0 -10.0 40]"; index = 1; }; / * Note that obj f i l e s have a .obj ending and s t r i p f i l e s have a .objf ending also, i f s t r i p f i l e s of the form <stripfileroot><level>.objf don't exist then a <objf i l e r o o t x l e v e l > . o b j w i l l be created so that ' s t r i p e ' can be run. */ /* Object bowl Loop { scale = [ 0.05 0.05 0.05 ]; of f s e t = [ 0 0 0 ] ; o b j f i l e r o o t = "data/flatbowl"; s t r i p f i l e r o o t = "data/flatbowl"; l e v e l s = 4; f i x e d = true; texture = "data/woodbowl.jpg"; t e x t u r e T i f f F i l e = "data/woodbowl.tif"; txtFileName = "data/woodbowl.txt"; }; */ Object rattleback Bezier { scale = [ 0.6 0.2 0.2 ]; rotate = -0.174 [ 0 0 1 ]; of f s e t = [0 0.0 0.2]; mass = 5; i n e r t i a = [ 1 0 0 ] [ 0 10 f i x e d = f a l s e ; controlp = [ -1.5 -1.5 -2.0] [ -0.5 -1.5 [ -1.5 -0.5 -1.0] [ -0.5 -0.5 [ -1.5 0.5 -1.0] [ -0.5 0.5 [ -1.5 1.5-2.0] [ -0.5 1.5 0 ] [ 0 0 10 ] ; -1.0] [ 0.5 0.0] [ 0 . 5 0.0] [ 0 . 5 -1.0] [ 0.5 -1.5 -0.5 0.5 1.5 -1.0] 0.0] 0.0] -1.0] [ 1.5 [ 1.5 [ 1.5 [ 1.5 -1.5 -2.0] -0.5 -1.0] 0.5 -1.0] 1.5 -2.0] 88 Object floor Bezier { scale = [ 1.0 1.0 offset = [0 0 0]; fixed = true; controlp = [ -1.5 -1.5 0.0] [ -0 [ -1.5 -0.5 0.0] [ -0 [ -1.5 0.5 0.0] [ -0 [ -1.5 1.5 0.0] [ -0 Contact C2 { objectA = floor; objectB = rattleback; codeA = "oponoponoponoponoponopono"; codeB = "onnononno"; / / tap start / / coords = 0.15 0.5 0.5 0.5 0; / / dcoords = 0 0 0 0 0 ; / / spin start coords = 0.45 0.45 0.5 0.5 0; dcoords = 0 0 0 0 - 1 ; / / this works... but be sure to use / / wmuxy = 0 and / / wmuz = small (0 or -0.0001 is goo) }; B.2 Subdiv is ion Surface F i l e Subdivision surface files use the obj file format. Here is the definition of the bowl seen on the left in Figure 4.3. v -1.0 0.0 1.0 v 0.5 -0.866 1.0 v 0.5 0.866 1.0 v -1.0 0.0 0.0 v 0.5 -0.866 0.0 v 0.5 0.866 0.0 v -2.0 0.0 -1.0 v 1.0 -1.732 -1.0 v 1.0 1.732 -1.0 1-0 ] ; .5 -1.5 0.0] .5 -0.5 0.0] .5 0.5 0.0] .5 1.5 0.0] [ 0.5 [ 0.5 [ 0.5 [ 0.5 -1.5 0.0] -0.5 0.0] 0.5 0.0] 1.5 0.0] [ 1.5 [ 1.5 [ 1.5 [ 1.5 -1.5 0.0] -0.5 0.0] 0.5 0.0] 1.5 0.0]; 89 f 9 8 7 f 4 5 6 f 4 6 3 f 3 1 4 f 6 5 2 f 2 3 6 f 5 4 1 f 1 2 5 f 2 1 7 f 7 8 2 f 1 3 9 f 9 7 1 f 3 2 8 f 8 9 3 90
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Fast contact evolution for piecewise smooth surfaces
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Fast contact evolution for piecewise smooth surfaces Kry, Paul G. 2000
pdf
Page Metadata
Item Metadata
Title | Fast contact evolution for piecewise smooth surfaces |
Creator |
Kry, Paul G. |
Date Issued | 2000 |
Description | Dynamics simulation of smooth bodies in contact is a critical problem in physically based animation and interactive virtual environments. We describe a technique which uses reduced coordinates to evolve a single continuous contact between Loop subdivision surfaces. The incorporation of both slip and no-slip friction into our algorithm is straightforward. The dynamics equations, though slightly more complex due to the reduced coordinate formulation, can be integrated easily using explicit integrators without the need for constraint stabilization. The use of reduced coordinates also confines integration errors to lie within the constraint manifold which is preferable for visualization. Our algorithm is suitable for piecewise parametric or parameterizable surfaces with polygonal domain boundaries. Because a contact will not always remain in the same patch, we demonstrate how a contact can be evolved across patch boundaries. We also address the issue of non-regular parameterizations occurring in Loop subdivision surfaces through surface replacement with n sided S-patch surfaces. Three simulations show our results. We partially verify our technique first with a frictionless system and then with a blob sliding and rolling inside a bowl. Our third simulation shows that our formulation correctly predicts the spin reversal of a rattleback top. We also present timings of the various components of the algorithm. |
Extent | 3905997 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-13 |
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.0051282 |
URI | http://hdl.handle.net/2429/10690 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-0450.pdf [ 3.73MB ]
- Metadata
- JSON: 831-1.0051282.json
- JSON-LD: 831-1.0051282-ld.json
- RDF/XML (Pretty): 831-1.0051282-rdf.xml
- RDF/JSON: 831-1.0051282-rdf.json
- Turtle: 831-1.0051282-turtle.txt
- N-Triples: 831-1.0051282-rdf-ntriples.txt
- Original Record: 831-1.0051282-source.json
- Full Text
- 831-1.0051282-fulltext.txt
- Citation
- 831-1.0051282.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-0051282/manifest