Controllable, Non-Oscillatory Damping for Deformable Objects by Herbert David Young B.Sc., The University of British Columbia, 2005 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in The Faculty of Graduate Studies (Computer Science) The University of British Columbia December 2007 c© Herbert David Young, 2007 Abstract This thesis presents a new method for the controllable damping of de- formable objects. The method evolves from physically based techniques; however, it allows for non-physical, but visually plausible motion. This flex- ibility leads to a simple interface, with intuitive control over the behaviour of the material. This method is particularly suited for strongly damped materials, which account for the majority of objects of interest to animation, since it produces non-oscillatory behaviour. This is similar to critical damping, except that it affects all modes independently. The new method is based on the minimiza- tion of a slightly modified version of total energy. This framework can be used to simulate many other physical phenomena, and therefore lends itself to coupling with other simulations. Implementation details for a simple example are given. Results are shown for varying parameters and compared to those produced by a traditional method. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Elasticity . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Energy Minimization . . . . . . . . . . . . . . . . . . 5 1.1.3 Perception . . . . . . . . . . . . . . . . . . . . . . . . 8 1.2 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 iii Table of Contents 2.1 Total Energy Minimization . . . . . . . . . . . . . . . . . . . 10 2.1.1 Modified Total Energy . . . . . . . . . . . . . . . . . 13 2.1.2 Explanation of τ . . . . . . . . . . . . . . . . . . . . . 14 2.2 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.1 Critical Damping of All Modes . . . . . . . . . . . . . 17 2.2.2 Equivalent Partial Differential Equation . . . . . . . . 18 2.2.3 New Properties . . . . . . . . . . . . . . . . . . . . . 21 3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1 Input and Pre-computation . . . . . . . . . . . . . . . . . . . 23 3.2 Each Time Step . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2.1 Minimization . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.2 Calculation of Total Energy . . . . . . . . . . . . . . 26 3.2.3 Update . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Pseudocode . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.1 New Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.1.1 Distance from the Rest State . . . . . . . . . . . . . . 32 4.1.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2 Traditional Method . . . . . . . . . . . . . . . . . . . . . . . 35 4.2.1 Choosing Damping Parameters . . . . . . . . . . . . . 36 4.2.2 Material Limits to Critical Damping . . . . . . . . . . 40 4.3 Comparison of Methods . . . . . . . . . . . . . . . . . . . . . 41 iv Table of Contents 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.1.1 Oscillation Control . . . . . . . . . . . . . . . . . . . 46 5.1.2 Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.1.3 Extending the Energy Minimization Framework . . . 47 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 v List of Tables 4.1 Parameters for Examples. . . . . . . . . . . . . . . . . . . . . 35 (a) Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . 35 (b) Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2 Damping Parameters for Traditional Method. . . . . . . . . . 42 (a) Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . 42 (b) Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . 42 vi List of Figures 2.1 Examples of under-, over-, and critical damping. . . . . . . . 16 2.2 The separate modes of two over-damped systems. . . . . . . . 16 (a) c = 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 (b) c = 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.1 Example 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 (a) Rest configuration . . . . . . . . . . . . . . . . . . . . . 33 (b) Initial displacement with fixed nodes marked with red dots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Example 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 (a) Rest configuration . . . . . . . . . . . . . . . . . . . . . 34 (b) Initial displacement with fixed nodes marked with red dots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.3 Method Comparison for Example 1. . . . . . . . . . . . . . . 43 (a) New Method. . . . . . . . . . . . . . . . . . . . . . . . . 43 (b) Traditional Method. . . . . . . . . . . . . . . . . . . . . 43 4.4 Method Comparison for Example 2. . . . . . . . . . . . . . . 44 (a) New Method. . . . . . . . . . . . . . . . . . . . . . . . . 44 (b) Traditional Method. . . . . . . . . . . . . . . . . . . . . 44 vii List of Programs 3.1 Pseudocode summary of algorithm . . . . . . . . . . . . . . . 30 viii Acknowledgements First and foremost, I thank my supervisor, Robert Bridson. His help, en- couragement, and inspiration have made this thesis possible. I am also thankful to Michiel van de Panne for his helpful comments, and fellow Imager students Christopher, Mike, and Stelian for their suggestions and motivation. As well, financial support from the Natural Sciences and Engineering Research Council of Canada allowed me to pursue this work. Finally, I would like to thank my family for their constant encouragement and support. For that, I am eternally grateful. ix To my wife Alison, for everything. x Chapter 1 Introduction Modelling elasticity is required for many simulations or animations. This is because elastic materials are both common and interesting. Most of the everyday objects we interact with have elastic properties, though some are more obvious than others. Many of the materials of interest to animation are elastic and strongly damped, such as hair, flesh, and clothing. With the growing popularity of computer-animated feature films and special effects, computer graphics are being pushed ever further. Initially, computer-generated animations were quite rigid, as deformable objects re- quired either a significant effort by the animators or complex computer sim- ulations. Since then, simulation techniques have become more developed and efficient, and with the ever advancing power of hardware, the capabili- ties of computer-generated content is always increasing. In order to be more realistic, or even more artistic, computer-generated content needs to model elastic objects. Simulations for film have very different requirements than those for engi- neering or research. When the target is entertainment, emphasis is placed on visual quality, speed of simulation, and ease of use, while scientific purposes demand proven accuracy. 1 Chapter 1. Introduction In this research, a method is presented which allows for intuitive control over the damping of an elastic material. The new method is inspired by physics, but motivated to offer a simple interface to the artist. The results, we believe, are physically plausible, even when the user requests physically impossible motion. 1.1 Related Work Several different topics in computer graphics provide important ground work for this new approach. Some cover previous methods which accomplish similar tasks, others provide motivation, and a few explain the choices made while developing this method. 1.1.1 Elasticity Deformable objects in general have been studied in computer graphics for over three decades, and physically based methods to simulate such objects have been in the field for at least two decades. For an overview of these methods and their history, the complementary survey papers of Gibson and Mirtich [7] and Nealen et al. [13] are excellent resources. Mass-Spring Systems One of the earliest of the physically based models was the mass-spring sys- tem. An array of point-masses — often regularly spaced in one, two, or three dimensions — are connected by simple, often linear, springs. Note that multiple connections of springs are often needed to keep the object from 2 Chapter 1. Introduction shearing, and that these springs usually include dampers as well. With a sufficient number of point-masses, the system can simulate complex shapes and deformations. These systems have been used in the early work of facial animation [21, 26], and have been used extensively for cloth [2, 3, 5]. The mass-spring system can also be used in other interesting ways: by dynami- cally adjusting the rest length of different springs, the model reacts to the changes as if being controlled by flexing muscles [12, 25]. Mass-spring models are intuitive, simple to implement, and fast to sim- ulate; however, the exact simulation of an elastic material is impossible [6]. The results depend on the number and distribution of point-masses, and special care must be taken when choosing individual spring parameters in order to even approximate a desired behaviour [6, 13]. Finite Element Methods A popular method in several fields, including computer graphics, for mod- elling elasticity is the finite element method (FEM). In general, it is used to solve partial differential equations. In terms of elastic simulation, the differential equation comes from continuum mechanics, ρẍ = ∇ · σ + f (1.1) where ρ is density, ẍ is the second time-derivative of position x, σ is stress, and f is external force [13]. FEM works by approximating the differential equation with a series of elements, often triangles or quadrilaterals in two dimensions. The approximation of each element is determined by the values 3 Chapter 1. Introduction associated with the nearby nodes and the basis functions. The finite element method has been used to simulate elasticity in com- bination with more complicated phenomena, such as brittle or ductile frac- ture [16, 17]. Even Green strain, with non-linear elastic forces, can be computed relatively easily with this method, which allows for greater defor- mation than linear elasticity. The accuracy depends on the choice of bases as well as the number of elements. The finite element method, which considers the object as a continuum, is more accurate than the mass-spring approach, which uses a discrete object model [7]. Unlike the mass-spring model, as the resolution is increased, and therefore the size of the elements is decreased, the model converges on a solution. Damping Parameters Both of the above methods still require damping parameters which do not easily relate to the time the system takes to come to rest. For the mass-spring method, choosing physically correct damping parameters is very difficult, if not impossible, since the discretization of the object affects the results of the damping parameters. Simply choosing parameters which critically damp each spring does not result in the object being critically damped. Even for other methods which can approximately use physical damping parameters, the correct values can be hard to find for certain materials, and the process of determining the correct values to obtain a certain result is mostly trial- and-error [7]. Modal dynamics have a damping coefficient for each independent mode; 4 Chapter 1. Introduction however, this is only efficient with a small number of modes of deformation, and is not extendable to plasticity [20]. Quasi-Static Methods Quasi-static motion ignores inertia and assumes that the material instanta- neously moves into a state of equilibrium after any change to the external forces. The method solves for the internal forces which will exactly match the external forces and keep the material in equilibrium. This method works well for materials where such an assumption about the inertia is approximately true, such as flesh in most circumstances [8, 24]. It allows simulations to run quickly, since the differential equation used by these methods no longer in- volves time, and a large amount of information can be pre-computed. More- over, estimating and tuning damping parameters is unnecessary, as damping and oscillation behaviour has been removed from the dynamics. However, the quasi-static assumption means that these methods are not suitable for ballistic motion, transient waves, creep, or other such phenom- ena, where inertia does produce visible results. This is because there is no temporal coherence; each step of the algorithm produces a result which is in equilibrium, ignoring the dynamics and the time the system would actu- ally take to reach this equilibrium point. The biggest limitation is that this assumption is not valid for many of the elastic objects commonly animated. 1.1.2 Energy Minimization Minimizing different energy terms can be a powerful method for physical simulation, and is often closely linked to how physicists understand these 5 Chapter 1. Introduction systems. Optimization can be used to find only certain forces, or it can solve for the entire system’s new configuration. Minimizing a sum of energies, ranging from gravitational and elastic po- tential to an artificial repulsion energy, can be used to enforce inter-particle constraints in cloth simulation [2]. It can also be used when creating pat- terns for apparel, where stretching and shearing need to be minimized [18]. In Lagrangian dynamics, one can minimize the integral of kinetic energy minus potential energy; however, this applies only to constrained systems, whereas we are interested in unconstrained dynamics. Quasi-static methods are essentially a minimization of the potential en- ergy, since the minimum of potential must be an equilibrium state. Rigid body contact can be phrased as a kinetic energy minimization with con- straints [11]. The pressure projection step used in fluid simulation can also be completed by minimizing kinetic energy, and because of the common framework, it has been used to couple fluids and solids [1]. Optimization of an energy-like term is also possible. For example, the minimization of an energy-like function has led to a new time integration scheme for dynamic systems [10]. Optimal control is another example, where a measure of control energy and the error from a target configuration is optimized to derive a control policy [4]. Damping Another physical phenomenon which can be treated with the minimization of energy framework is damping. By adding a regularization term to kinetic energy, the optimization method is equivalent to using a Backward Euler 6 Chapter 1. Introduction integration scheme. To see this, start with a typical FEM discretization for force, F = DS (1.2) where F is a vector with component Fi being the net force on modal vari- able i, S is the stress tensor at each quadrature point, and D is a discrete approximation to the divergence operator, mapping stress at the quadrature points to force on the nodes. If we consider just the damping stress, then S = −ADTv (1.3) where A is a symmetric positive definite matrix which encodes the damping parameters, v is the velocity, and DTv is an approximation of strain rate. A Backward Euler solve for damping is Mvnew = Mṽ +∆tDS (1.4) S = −ADTvnew (1.5) with M being the mass matrix and ṽ is the intermediate velocity with non- damping forces already integrated. The usual method to solve this is to substitute (1.5) into (1.4), vnew = ( M+∆tDADT )−1 Mṽ (1.6) 7 Chapter 1. Introduction but if we minimize kinetic energy with a regularization term, ∆t2 S TA−1S: min S 1 2 vTnewMvnew + ∆t 2 STA−1S (1.7) ⇔ 0 = ( dvnew dS )T Mvnew +∆tA−1S (1.8) 0 = ∆tDTM−1Mvnew +∆tA−1S (1.9) S = −ADTvnew (1.10) then we get the same result as a Backward Euler step. 1.1.3 Perception A very important area of discussion in computer graphics is determining how much accuracy is truly necessary. Physically correct motion is important in order to create believable animations. However, if the user is unable to see the difference between a perfectly accurate simulation and a good approximation, then the latter may be more desirable. A method which produces an acceptable approximation may have several benefits, such as being easier to use, faster to implement, or more efficient to run. Studies have shown people’s ability to understand physical systems is limited in various ways. Despite the deterministic behaviour of most physi- cal systems, if more than one parameter or dimension is significant, people lose their intuition [9, 22]. Perceptual studies have been done on deter- mining exactly what motion appears plausible to the viewer, but a robust and accurate metric is still unavailable [19]. A recent study suggests that predicting the behaviour of a physical system and determining plausibility 8 Chapter 1. Introduction are independent cognitive processes, and that approximations can be made without producing visually objectionable results [15]. 1.2 Thesis Overview First, the new method is introduced in Chapter 2, and the major compo- nents are explained. In Chapter 3, an implementation of the method is described, and details specific to implementing are provided. The results of the implementation are then discussed and demonstrated in Chapter 4. Finally, Chapter 5 concludes with what has been accomplished, and what further work this could lead to. 9 Chapter 2 Method Optimization has been used to simulate many different physical systems, as discussed in Section 1.1.2. In some circumstances — including damp- ing with a regularization term, rigid-body contact, and incompressible fluid flow — dynamics can be modelled by minimizing kinetic energy, and quasi- static systems by minimizing potential energy. In order to capture both the dynamic behaviour of the former model and the control-like treatment of damping in the latter, we start by combining the concepts in a mini- mization of total energy. The resulting method is dependant on the time discretization. To overcome this problem, we use a modified version of total energy. 2.1 Total Energy Minimization If we simply minimize the total energy, the sum of kinetic and potential en- ergies, with respect to stress, or force, then an undesirable property results: the method’s apparent damping rate is reduced as the step size becomes smaller. To demonstrate this, we consider an elastic bar in one-spatial dimension. After taking a Backward Euler step, we use Fourier analysis and minimize 10 Chapter 2. Method the total energy with respect to force. First, define transforms for displace- ment un, velocity vn, and force f , at time step n, un = Ane √−1kx (2.1) vn = Bne √−1kx (2.2) f = Fe √−1kx (2.3) where x is the material space, or rest, position and k is the wave number. The time discretization is then Bn+1 = Bn + ∆t ρ F (2.4) An+1 = An +∆tBn+1 (2.5) = An +∆tBn + ∆t2 ρ F (2.6) Total energy is the sum of kinetic and potential energy over the domain of the object, Ω, which in one dimension is the object length. The energy terms will be defined as KEn = ∫ Ω 1 2 ρv2n (2.7) PEn = ∫ Ω 1 2 E ( ∂un ∂x )2 (2.8) TEn = KEn + PEn (2.9) where ρ is the material density, E is Young’s modulus of elasticity, and ∂u∂x is the strain. The approximate effect on the displacements follows from minimizing 11 Chapter 2. Method the total energy at the next time step, min ∫ Ω KEn+1 + PEn+1 (2.10) ⇔ min ∫ Ω 1 2 ρv2n+1 + 1 2 E ( ∂un+1 ∂x )2 (2.11) ⇔ min F ∫ Ω ρ 2 ( Bn + ∆t ρ F )2 ∣∣∣e√−1kx∣∣∣2 + ∫ Ω E 2 ( An +∆tBn + ∆t2 ρ F )2 ∣∣∣√−1ke√−1kx∣∣∣2 (2.12) ⇔ min F Ωρ 2 ( Bn + ∆t ρ F )2 + |k|2ΩE 2 ( An +∆tBn + ∆t2 ρ F )2 (2.13) Finding the minimum is equivalent to taking the derivative of this objective function with respect to F , and finding the zero, Ω∆t ( Bn + ∆t ρ F ) + |k|2ΩE∆t 2 ρ ( An +∆tBn + ∆t2 ρ F ) = 0 (2.14) and therefore F is ∆t ρ ( 1 + |k|2E∆t 2 ρ ) F = −|k|2E∆t ρ An − ( 1 + |k|2E∆t 2 ρ ) Bn (2.15) F = − |k| 2E 1 + |k|2E∆t2ρ An − ρ∆tBn (2.16) Applying this to (2.4) and (2.6) gives Bn+1 = − |k| 2E∆t ρ+ |k|2E∆t2An (2.17) An+1 = 1 1 + |k|2E∆t2ρ An (2.18) 12 Chapter 2. Method and therefore, with the total time T = n∆t, An = ( 1 1 + |k|2E∆t2ρ )n A0 (2.19) An ≈ e−n|k| 2 E∆t2 ρ A0 (2.20) A(T ) ≈ e−|k|2 Eρ T∆tA0 (2.21) Since the behaviour depends on the step size, the method becomes less efficient as ∆t → 0. This is an undesirable property for several reasons. Smaller values of ∆t mean smaller steps in time, and therefore more steps are needed than with a larger ∆t to simulate the same total time. If the method becomes less efficient, then the time to simulate will grow faster than the number of steps needed, making very small values of ∆t impracticable. 2.1.1 Modified Total Energy In order to avoid this problem, we minimize a slightly different value. We multiply the potential energy by a dimensionless fraction, which means we maintain the units of energy. This new coefficient has the units of time over time which results in not only eliminating the above problem, but also giving the user control over the length of time the simulation takes to come to rest. Our new minimization has the form min ∫ KE + τ ∆t PE (2.22) Now by adjusting τ , we have control over the time the system takes to come to rest, which we will call trest. The effect of τ is inversely related to 13 Chapter 2. Method trest, which will become more clear with a similar derivation to (2.21), but with the new coefficient: min ∫ Ω 1 2 ρ(vn+1)2 + τ 2∆t E ( ∂un+1 ∂x )2 (2.23) F = − |k| 2E τ∆t 1 + τ |k|2E∆tρ An − ρ∆tBn (2.24) An = ( 1 1 + τ |k|2E∆tρ )n A0 (2.25) An ≈ e−nτ |k| 2 E∆t ρ A0 (2.26) A(T ) ≈ e−τ |k|2 Eρ TA0 (2.27) Therefore, if the system comes to rest by time T = trest1 using τ = τ1, then by using τ = ατ1, the same rest state is achieved by time 1α t rest 1 . This gives the user a much simpler interface to control the behaviour of elastic material than manually adjusting damping parameters. It is important to note that the minimization of the modified total energy does not have the same inefficiency problem as the unmodified version. In (2.27), unlike (2.21), there is no ∆t, and therefore the size of the time step does not affect the rate at which the system comes to rest. 2.1.2 Explanation of τ The behaviour of τ seems unintuitive at first, especially considering it has the units of time and not the inverse. Initially, the coefficient of potential energy, τ ∆t , looks like the number of steps of size ∆t to reach time τ . However, as we increase this fraction, the time to reach rest, or any particular state, 14 Chapter 2. Method decreases. So using the same ∆t and increasing τ , decreases the time to come to rest. To understand τ better, we will take a look at the damping of a simple mass-spring system. Taking the basic equation of motion for displacement x(t), a function of time, with mass m, damping coefficient c, and stiffness k, we get mẍ+ cẋ+ kx = 0 (2.28) where ẋ = dxdt and ẍ = d2x dt2 . Assuming that the solution to this ordinary differential equation is of the form x(t) = ert, we can find r by solving mr2 + cr + k = 0 (2.29) using the quadratic formula, r = −c±√c2 − 4mk 2m (2.30) In this form, (2.30), we can easily identify how the system is damped by the value in the square root, α = c2 − 4mk. If α is zero, so c2 = 4mk, then the system is critically damped, and there exists one, real solution: r1 = r2 ∈ R. In the over-damped case, α is greater than zero and there are two unique, real solutions: r1 6= r2 ∈ R; in the under-damped case, α is less than zero and there are two unique, complex solutions: r1 6= r2 ∈ C. By graphing an example for each case, Figure 2.1, the effects of different values of α are very clear. To solve (2.28), we use x(t) = Aer1t+Ber2t. These graphs use m = 1, k = 1, A = 0.5, B = 0.5 with c = 1 for under-damping, 15 Chapter 2. Method 0 2 4 6 8 10 12 14 16 −0.5 0 0.5 1 rest under−damped critically damped over−damped Figure 2.1: Examples of under-, over-, and critical damping. c = 2 for critical damping, and c = 4 for over-damping. The values of r in each case are −0.5 ± √ 3 2 i for under-damped, −1 for critically damped, and −2±√3 for over-damped. 0 2 4 6 8 10 12 14 16 −0.2 0 0.2 0.4 0.6 0.8 1 (a) c = 3 0 2 4 6 8 10 12 14 16 −0.2 0 0.2 0.4 0.6 0.8 1 (b) c = 6 Figure 2.2: The separate modes of two over-damped systems. Since the differential equation (2.28) is second-order, it involves finding the solution of a quadratic equation. The two roots mean that there are two modes in the solution. In the critically damped case, the two modes are equal, but in the other cases, one mode comes to rest faster; let that 16 Chapter 2. Method time be denoted as tfast and the other slower time as tslow. If the damping parameter is increased past critically damping, then tslow increases and so does the time the whole system takes to come to rest; however, as tslow increases, tfast decreases. This can be observed in Figure 2.2 where the modes of an over-damped system are plotted separately. As shown in Section 2.2.1, our new method of simulating elastic mate- rial critically damps all modes, and because of this, the duality of the above method does not occur. The minimization works on each mode indepen- dently, and since the minimum will be related to tfast, the other mode, tslow is not seen. However, if tslow is theoretically increased, the system comes to rest sooner, and it is this reason, that τ can be treated as the theoretical value of tslow, with units of time and affecting the material inversely. The behaviour of τ can also be further understood through the analysis given in Section 2.2.2. 2.2 Analysis 2.2.1 Critical Damping of All Modes A major benefit of this method is the damping control. As mentioned above, the user can adjust the time the system takes to come to rest with a single, easily-understood parameter. Another effect of the minimization framework is that all modes are critically damped. Such behaviour is not attainable with traditional methods, regardless of the choice of damping parameters. This unique result of the new method is due to the behaviour of mini- 17 Chapter 2. Method mization. Parseval’s theorem states that ∫ Ω |A|2 = c(Ω) ∞∑ n=−∞ |an|2 (2.31) where an are the coefficients in the Fourier series of A and c is some constant depending on Ω. If we apply this to our modified total energy definition in (2.23), we get min σ ρ 2 ∫ Ω (vn+1)2 + τE 2∆t ∫ Ω ( ∂un+1 ∂x )2 (2.32) minbσ ρ 2 c(Ω) ∞∑ k=−∞ |v̂n+1k |2 + τE 2∆t c(Ω) ∞∑ k=−∞ ∣∣ûn+1k √−1k∣∣2 (2.33) where v̂n+1, ûn+1, and σ̂ are the Fourier series coefficients for vn+1, un+1, and σ. Since minimizing with respect to the vector σ̂ is equivalent to minimizing with respect to each component σ̂k independently, each mode of energy depends on only one mode of the stress, and therefore each mode is critically damped. 2.2.2 Equivalent Partial Differential Equation To better understand the new method, we derive the partial differential equation that it is essentially solving. To simplify the process, we show only the one-dimensional problem with spatially uniform material properties — with displacement u, velocity v, stress σ, and strain ε = ∂u∂x — and we do not expand strain, even though it is not a primary variable. The problem is 18 Chapter 2. Method discretized in time using vn+1 = vn + ∆t ρ ∂σ ∂x (2.34) un+1 = un +∆tvn+1 (2.35) and therefore εn+1 = εn +∆t ∂vn+1 ∂x (2.36) We now use variational calculus and assume that we have a stress, σ, which minimizes the modified total energy, given in (2.22). If we perturb that solution by δs, for some arbitrary stress s, such that the stress of the new system is σ + δs, then the new energy can be written as a function of δ: g(δ) = ∫ ρ 2 ∣∣∣∣vn+1 + ∆tρ δ ∂s∂x ∣∣∣∣2 + τ∆t ∫ E 2 ∣∣∣∣εn+1 + ∆t2ρ δ ∂2s∂x2 ∣∣∣∣2 (2.37) Since we assume that σ minimizes the energy, δ = 0 must be the minimum for the new function, g(δ). At this minimum, the derivative of the function must be zero: g′(δ) = ∫ ρ 2 ( 2 ∆t ρ vn+1 ∂s ∂x + 2δ ( ∆t ρ ∂s ∂x )2) + τ ∆t ∫ E 2 ( 2 ∆t2 ρ εn+1 ∂2s ∂x2 + 2δ ( ∆t2 ρ ∂2s ∂x2 )2) (2.38) g′(0) = ∫ ∆tvn+1 ∂s ∂x + τ ∫ E ∆t ρ εn+1 ∂2s ∂x2 = 0 (2.39) Considering the domain between a and b, this last equation can be sim- 19 Chapter 2. Method plified by dividing by ∆t and using integration by parts, 0 = ∫ b a vn+1 ∂s ∂x + τE ρ ∫ b a εn+1 ∂2s ∂x2 (2.40) = − ∫ b a ∂vn+1 ∂x s+ [ vn+1s ]b a + τE ρ (∫ b a ∂2εn+1 ∂x2 s+ [ εn+1 ∂s ∂x ]b a − [ ∂εn+1 ∂x s ]b a ) (2.41) Ignoring boundary conditions to further simplify the equation results in ∫ ( −∂v n+1 ∂x + τE ρ ∂2εn+1 ∂x2 ) s = 0 (2.42) Since s is an arbitrary stress, the only way the above equation holds true is when the first factor of the integrand is exactly zero, −∂v n+1 ∂x + τE ρ ∂2εn+1 ∂x2 = 0 (2.43) Using (2.36), this is equivalent to εn+1 − εn ∆t = τE ρ ∂2εn+1 ∂x2 (2.44) which is Backward Euler applied to a strain diffusion equation, ∂ε ∂t = τE ρ ∂2ε ∂x2 (2.45) Using this new equation in Fourier space shows the same results as the previous section: each mode is independently and critically damped. With 20 Chapter 2. Method ε̂i being the weight for mode i of strain, dε̂i dt = −τE ρ i2ε̂i (2.46) and therefore the decay of each mode is exponential, ε̂i(t) = e − τE ρ i2t ε̂i(0) (2.47) Since strain is the spatial derivative of displacement, taking the antideriva- tive with respect to i leads to similar equations for displacement. A similar analysis in two dimensions leads to two diffusion equations for the displacement, 1 τ ∂u1 ∂t = ( λ+ 2µ ρ ) ∂2u1 ∂x2 + ( λ+ 2µ ρ ) ∂2u2 ∂x∂y + ( 2µ ρ ) ∂2u1 ∂y2 (2.48) 1 τ ∂u2 ∂t = ( 2µ ρ ) ∂2u2 ∂x2 + ( λ+ 2µ ρ ) ∂2u1 ∂x∂y + ( λ+ 2µ ρ ) ∂2u2 ∂y2 (2.49) both analogous to Equation (2.45). Note that λ and µ, the Lamé parameters, are both functions of E, Young’s modulus of elasticity, and ν, Poisson’s ratio. They 2.2.3 New Properties Choosing the stress by optimization leads to a new constitutive law. The relationship between stress and strain is no longer defined physically, which gives the method many of its characteristics. 21 Chapter 2. Method Regular elasticity leads to a wave equation with a damping term, D, ∂2ε ∂t2 = E ρ ∂2ε ∂x2 +D (2.50) but the new method solves a diffusion equation, (2.45), which can be written as 1 τ ∂ε ∂t = E ρ ∂2ε ∂x2 (2.51) Therefore the strain energy diffuses over the whole object without prop- agating in waves, which is why the new method shows critical damping behaviour. The new dynamics are non-oscillatory and the rate of strain diffusion can be controlled by the new parameter τ . The higher the value of τ , the faster strain diffuses, and therefore the faster the system comes to rest. By taking the previous energy-derived techniques and finding a middle ground, this new method has several interesting properties. The kinetic en- ergy term captures inertial effects, and since it still solves equation (1.1), momentum is conserved and ballistic motion is possible. The augmented potential energy term takes the instant damping of quasi-static motion and makes it a controllable behaviour. The system can be more directly con- trolled than previously possible, and the results are plausible even if the requested action is physically impossible. 22 Chapter 3 Implementation This section provides the details of a two-dimensional, triangle-based imple- mentation. This simple implementation assumes the simulation involves a single, homogenous, isotropic material. No collision detection is involved, including that for self-intersection. 3.1 Input and Pre-computation The inputs are a triangle mesh, along with a list of fixed vertices and a list of initial vertex displacements, and the global material parameters of the ob- ject. These include the mass associated with each vertex, Lamé parameters λ and µ, as well as the value of the new parameter τ . The deformable object, defined in its original rest state using material coordinates p, has world space coordinates of x. The strain on each element, a triangle in this case, is calculated with Green strain, εij = 1 2 ( ∂ui ∂xj + ∂uj ∂xi + ∂ui ∂xj ∂uj ∂xi ) (3.1) where ui is the displacement at vertex i. This strain tensor can also be 23 Chapter 3. Implementation written as G = 1 2 (ATA− I) (3.2) where A = ∂x∂p is the deformation gradient and I is the identity matrix. As shown by Teran et al. [23], if we define the columns of two, 2 × 2 matrices, Dm and Ds, using the edge vectors of each triangle in both material and world space coordinates, Dm = ( p1 − p0 ∣∣ p2 − p0 ) (3.3) Ds = ( x1 − x0 ∣∣ x2 − x0 ) (3.4) then Ds = ADm (3.5) A = DsD−1m (3.6) The advantage with this definition is that theD−1m matrix is constant and can be pre-computed, and since Ds is simple to construct, the computation of the strain, using equation (3.2) with (3.6), is fast and easy. 3.2 Each Time Step For each time step, we find the stress tensors which minimize the modified total energy. Then we calculate the force produced by that stress, and use these to update the position of each vertex. 24 Chapter 3. Implementation 3.2.1 Minimization Virtually any optimization algorithm can be used for this method. For this example we used steepest descent, because of the ease of implementation and to keep the main focus of this new method clear. For each iteration, steepest descent, or gradient descent, takes a step pro- portional to the negative gradient. For this method, the update at iteration i is σi+1 = σi − αi∇TE(σi) (3.7) and TE(σi+i) ≤ TE(σi) (3.8) where TE(σ) is the total energy at stress σ. Step Length The choice of αi is arbitrary, as long as the inequality (3.8) is satisfied for all i ≥ 0. Several line search strategies exist, each with a balance between being fast to compute and returning better results; see Nocedal and Wright for some examples [14]. For this implementation, we start with a default value of α, and if con- dition (3.8) fails, then the step length is halved. This is repeated until a value of α which satisfies the non-increasing condition is found, or until α is sufficiently small. This procedure does work; however, at some points, many iterations are needed because of unnecessarily slow progress. This is the result of the default step length being too small. 25 Chapter 3. Implementation To reduce this problem, without turning to a more complex algorithm, the initial value of α is doubled whenever an iteration does not need to reduce the initial step size. If the previous iteration does reduce αi−1 from the initial guess, then the initial value is set to the final value of the previous iteration, αinitiali = αi−1. This adaptive method quickly moves over large areas with small gradients, but readily adapts to less smooth areas. It is similar to the technique used by Okabe et al. [18]. 3.2.2 Calculation of Total Energy The objective function used in the optimization is the modified definition of total energy, TE(σ) = KE(σ) + τ ∆t PE(σ) (3.9) and is calculated exactly as shown, with the kinetic and potential energy functions as defined below. Kinetic Energy The kinetic energy function first computes the force on each vertex resulting from the given stress. This force is used to update the velocity of each vertex, and then the kinetic energy of the whole system is given by KE(σ) = 1 2 v(σ)TMv(σ) (3.10) where M is the mass matrix, and velocity has been updated using vnew(σ) = vold +∆tM−1F(σ) (3.11) 26 Chapter 3. Implementation The old velocity is from the last time step, and each component of force, Fi, is the sum of forces on vertex i, from each adjacent triangle, Fi(σ) = ∑ t∈adj(i) σt (xtk − xtj )⊥ 2 (3.12) where adj(i) are all the triangles with vertex i, σt is the stress on triangle t, and xtk and xtj are the two vertices of triangle t opposite vertex i. Potential Energy The potential energy of the whole system is the sum of the potential energy of each triangle. The potential energy density on each triangle is constant, since the stress has been approximated as constant across each element. PE(σ) = ∑ t PEt(σt) (3.13) For a single triangle, the potential energy is PEt(σt) = ∫∫ Ωt 1 2 σt : G(σt) (3.14) where A : B = ∑ i,j AijBij is the double contraction operator. Combining our standard definition of stress, σij = λ 2∑ i=1 εiiδij + 2µεij (3.15) 27 Chapter 3. Implementation where δij is the Kronecker delta function defined as δij = 1 if i = j0 if i 6= j (3.16) with ε = G and (3.14), gives PEt(σt) = at 2 λ( 2∑ i=1 Gii )2 + 2µ 2∑ i=1 2∑ j=1 G2ij (3.17) The last line, (3.17), replaces the integral with a multiplication by the area of triangle t, at, since PEt(σt) is constant over each element. Gradient of Total Energy Now that we can compute the total energy at any given stress, we use finite differences to estimate the total energy gradient. Since stress is a symmetric 2 × 2 matrix, only three variables need to considered. As well, a change in the stress of one triangle only affects that triangle and its neighbours, making this a local operation and therefore feasible for a large simulation. Centred differences are used, resulting in the following approximation to the total energy gradient, ∇TE(σ) = 1 2δ TE(σ+∆11)−TE(σ−∆11) TE(σ+∆12)−TE(σ−∆12) TE(σ+∆21)−TE(σ−∆21) TE(σ+∆22)−TE(σ−∆22) (3.18) 28 Chapter 3. Implementation where δ is some small constant and ∆kl is a 2× 2 matrix defined as ∆klij = δ if ij = kl or ji = kl0 otherwise (3.19) 3.2.3 Update Once the stress that minimizes the objective function is found, it is used to update all, non-fixed nodes in the mesh. The the new velocities are calculated with (3.11), and the new positions by xnew(σ) = xold +∆tvnew(σ) (3.20) 3.3 Pseudocode The new method, as described above and as we implemented, is summarized in Program 3.1. 29 Chapter 3. Implementation Program 3.1 Pseudocode summary of algorithm // initialize read in mesh compute the inverse rest edge matrix for each triangle read in and apply initial displacement // steepest descent calculate gradTE with finite differences stress_trial = stress - alpha * gradTE calculate TE_new(stress_trial) while TE_new > TE_old alpha = alpha / 2 recalculate stress_trial and TE_new stress = stress_trial // update mesh for each non-fixed vertex calculate force(stress) v_new = v_old + dt/m * force x_new = x_old + dt * v_new 30 Chapter 4 Results 4.1 New Method The implementation described in Chapter 3 has been run with varying in- puts in order to observe the new method. In practice, the new approach to damping control is easy to use, and produces visually-plausible elastic motion. Ideally, the exact value of τ for a given time to come to rest would be calculated; however, realistically, τ is estimated, trest is observed, and then the desired value of τ is computed exactly. Although this leads to two simulations being needed, it is still a considerable improvement from the many iterations of trial-and-error commonly taken with traditional damping parameters. Theoretically, it would be possible to calculate τ beforehand, but because of the complex relationship between τ and trest, it is faster to run through the simulation an extra time. This action is further supported by the ability of the method to take large time steps, and to quickly deliver meaningful practice runs. 31 Chapter 4. Results 4.1.1 Distance from the Rest State In order to quantitatively measure when each simulation comes to rest, we have graphed two measures of the distance from the rest state. The first metric, d1, is the p-2, or Euclidean, norm of a matrix containing all the differences of each vertex’s position from the rest state, d1 = ‖x− p‖2 (4.1) and the second, d2 is the maximum distance any vertex is from its rest position, d2 = max i∈Ω ‖xi − pi‖2 (4.2) The first measure, d1, is global and considers how far the object as a whole is from rest, while d2 is local and measures the most extreme vertex. 4.1.2 Examples Example 1 The first example is a simple mesh with twelve triangles arranged in a trape- zoid, shown in its rest state in Figure 4.1a. The initial displacement, Figure 4.1b, is a single vertex at the top moved inward, and the vertices along the bottom are all fixed. The material parameters are given in Table 4.1a. Note that here, ρ is similar to the material density, but is just the mass at each vertex, and it is used to construct the mass matrix, M. In order to demonstrate the behaviour of the simulation, the distance from the rest state using the two metrics described in Section 4.1.1 are plot- 32 Chapter 4. Results ted against time. Figure 4.3a shows this example with two different values of τ . The graph demonstrates the effects of τ quite well: the simulation using τ = 1 comes to rest at approximately T = 3s, while τ = 4 arrives at the same state by T = 0.75s. Note that the behaviours of both simulations are the same qualitatively, and that they show a critically-damped-like approach to the rest state. This is because adjusting the value of τ is essentially a time-scaling operation. Also note that the two different metrics are virtually identical in this example because of the small number of vertices. (a) Rest configuration (b) Initial displacement with fixed nodes marked with red dots Figure 4.1: Example 1. Example 2 The second example is a mesh with 414 triangles laid out roughly in an annulus. The rest state of the mesh and the initial deformation is shown in Figure 4.2. In this example, the bottom vertices all fixed, and initially, the annulus is being pulled to the left and released. The material and time parameters for this animation are given in Table 4.1b. The values of ∆t are different for the two examples in order to simulate 33 Chapter 4. Results the larger mesh in a similar time, not because of any restrictions on size of the time steps. The behaviour over time, shown in Figure 4.4a, provides more interesting detail than the first example. The relationship between the rest times of τ = 1 and τ = 0.33 is as expected; however, τ = 2 does not quite come to rest in half the time of τ = 1. This is likely due to numerical error, as the rest times are becoming relatively small. Also of interest is the oscillatory behaviour of d2. The object as a whole still behaves as if critically damped, seen with d1, but individual vertices may still overshoot the rest state. This gives the actual animation an interesting motion, while still effectively damping the object as desired. (a) Rest configuration (b) Initial displacement with fixed nodes marked with red dots Figure 4.2: Example 2. 34 Chapter 4. Results Parameter Value λ 3 µ 1 ρ 1 ∆t 0.01s Tfinal 5s (a) Example 1 Parameter Value λ 2 µ 0.5 ρ 0.1 ∆t 0.05s Tfinal 5s (b) Example 2 Table 4.1: Parameters for Examples. 4.2 Traditional Method We have also implemented a two-dimensional, triangle-based, traditional technique for comparison. As much as possible, the same framework was used as described in Chapter 3. We start by computing the strain tensor with the same definition, (3.2), and same technique as we used with the new method. Then elastic stress, σe, is calculated by (3.15) and is used to find the elastic force, Fe, with (3.12). Intermediate velocities and positions are now given by ṽ = vold +∆tM−1Fe (4.3) x̃ = xold +∆tṽ (4.4) These intermediate values are used to calculate the strain rate, Ġ = 1 2 (ȦTA+AT Ȧ) (4.5) 35 Chapter 4. Results and Ȧ was calculated by Ȧ = ḊsD−1m , with Dm from (3.3) and Ḋs = ( ṽ1 − ṽ0 ∣∣ ṽ2 − ṽ0 ) (4.6) We can now find the damping stress, σd, similarly to (3.15), σdij = φ 2∑ i=1 Ġiiδij + 2ψĠij (4.7) and the damping force, Fd, is given by (3.12) using σd. The total force is the sum of the elastic and damping forces, F = Fe+Fd, and the final velocities and positions are given by vnew = vold +∆tM−1F (4.8) xnew = xold +∆tvnew (4.9) 4.2.1 Choosing Damping Parameters In order to reasonably compare the new method to existing techniques, we must be careful when selecting the damping parameters for the traditional methods. Since these methods can only critically damp a single mode, we find parameters which critically damp the mode which is most visually ob- vious, corresponding to either the size of the elements or the size of the object. 36 Chapter 4. Results The differential equation (1.1) without external force can be written as ρüi = ∂σij ∂xj (4.10) = ∂ ∂xj (λεkkδij + 2µεij − φε̇kkδij − 2ψε̇ij) (4.11) Define strain with εij = ( ∂ui ∂xj + ∂uj ∂xi ) (4.12) and the Fourier transform for displacement using u = U(t)e √−1kx (4.13) This is similar to (2.2), but continuous in time. Substituting (4.13) into (4.12) gives εij = Ui √−1kje √−1kx + Uj √−1kie √−1kx (4.14) and therefore ∂εij ∂xj = Ui( √−1kj)2e √−1kx + Uj( √−1ki)( √−1kj)e √−1kx (4.15) = −(Uik2j + Ujkikj)e √−1kx (4.16) 37 Chapter 4. Results Expanding this in three dimensions gives − (Uik2j + Ujkikj) = U1|k|2 + U1k21 + U2k1k2 + U3k1k3 U2|k|2 + U1k1k2 + U2k22 + U3k2k3 U3|k|2 + U1k1k3 + U2k2k3 + U3k23 (4.17) = |k|2 + k21 k1k2 k1k3 k1k2 |k|2 + k22 k2k3 k1k3 k2k3 |k|2 + k23 U1 U2 U3 = (|k|2I + k ⊗ k)U (4.18) Now the equation (4.11) can be written ρÜ+ ( µ|k|2I + (µ+ λ)(k ⊗ k))U+(ψ|k|2I + (ψ + φ)(k ⊗ k)) U̇ = 0 (4.19) The matrix |k|2I + k⊗ k has one eigenvector parallel to k and two more perpendicular to k. Let S be the component of U which is perpendicular to k, and T the component parallel to k. Now we can rewrite (4.19) with two independent equations, ρS̈ + µ|k|2S + ψ|k|2Ṡ = 0 (4.20) ρT̈ + (2µ+ λ)|k|2T + (2ψ + φ)|k|2Ṫ = 0 (4.21) These two simple differential equations, (4.20 & 4.21), can now be eas- ily solved. Assuming that the solution to (4.20) looks like S(t) = ert, we 38 Chapter 4. Results substitute this in and solve for r, 0 = ρr2 + ψ|k|2r + µ|k|2 (4.22) r = −ψ|k|2 ±√ψ2|k|4 − 4ρµ|k|2 2ρ (4.23) Since we want to find parameters which critically damp the material, we know that the discriminant of (4.23) must be zero, 0 = ψ2|k|4 − 4ρµ|k|2 (4.24) ψ = √ 4ρµ |k|2 (4.25) By doing the same with (4.21), we get ρr2 + (2ψ + φ)|k|2r + (2µ+ λ)|k|2 = 0 (4.26) and 0 = (2ψ + φ)2|k|4 − 4ρ(2µ+ λ)|k|2 (4.27) φ = √ 4ρ(2µ+ λ) |k|2 − 2ψ (4.28) Now, to critically damp any particular wavelength Λ, we choose the mode k = 2piΛ , and put this — along with material parameters λ, µ, and ρ — into equations (4.25) and (4.28). 39 Chapter 4. Results 4.2.2 Material Limits to Critical Damping It is interesting to note here what equations (4.25) and (4.28) imply. Ther- modynamics requires that the damping parameters be non-negative, so com- bining these equations and φ ≥ 0 means √ 4ρ(2µ+ λ) |k|2 − 2 √ 4ρµ |k|2 ≥ 0 (4.29) 2 √ ρ |k| (√ 2µ+ λ− 2√µ ) ≥ 0 (4.30)√ 2µ+ λ ≥ 2√µ (4.31) 2µ+ λ ≥ 4µ (4.32) λ ≥ 2µ (4.33) Both steps (4.31) and (4.32) are because ρ, |k|, µ, λ ≥ 0. This means that critical damping at a particular wavelength for both shear and compression waves is only possible in materials where the in- equality (4.33) holds. Substituting the definitions of λ and µ, µ = E 2(1 + ν) (4.34) λ = Eν (1 + ν)(1− 2ν) (4.35) into (4.33), and assuming that E > 0 and −1 < ν < 12 , this can be simplified 40 Chapter 4. Results further, Eν (1 + ν)(1− 2ν) ≥ 2 ( E 2(1 + ν) ) (4.36) ν (1− 2ν) ≥ 1 (4.37) ν ≥ 1 3 (4.38) Therefore, with assumptions which are true for all practical materials, critical damping is only possible when Poisson’s ratio, ν, satisfies (4.38). This is an interesting limit, as many metals have a ratio near one third, and several common materials have ratios less than a third, such as concrete, glass, and steel. 4.3 Comparison of Methods Now that we can choose damping parameters based on critically damping a desired wavelength, we will compare the new method to three different pairs of parameters. One set will critically damp the wavelength which corresponds to the size of the elements, the length of a triangle edge in this case, one will correspond to the radius, and the other to the diameter of the object. These parameters are given for both examples in Table 4.2. In order to compare the behaviours, the distance from rest is plotted for simulations using the new method with varying values of τ , and the traditional method for the above mentioned wavelengths. Figure 4.3 shows Example 1, where the traditional method shows much more oscillation than the new method, which is capable of critically damping both faster and 41 Chapter 4. Results slower than the other method. In the second example, Figure 4.4, the different wavelengths cause either too much oscillation overall or are over-damped. The new method shows the desired result for a heavily damped material, and although there are some oscillations shown in the d2 norm, both metrics show the new method coming to rest as fast or faster than the other method. Parameter Value Λ 1 k 2pi φ 0.053134 ψ 0.22508 Λ 2 k pi φ 0.10627 ψ 0.45016 Λ 4 k 12pi φ 0.21254 ψ 0.90032 (a) Example 1 Parameter Value Λ 1 k 2pi φ 0.022622 ψ 0.050329 Λ 8 k 14pi φ 0.18098 ψ 0.40263 Λ 16 k 18pi φ 0.36196 ψ 0.80527 (b) Example 2 Table 4.2: Damping Parameters for Traditional Method. 42 Chapter 4. Results 0 1 2 3 4 5 0 0.1 0.2 0.3 0.4 0.5 0.6 time (s) d 1 , d 2 τ = 1, d2 τ = 4, d2 τ = 1, d1 τ = 4, d1 (a) New Method. 0 1 2 3 4 5 0 0.1 0.2 0.3 0.4 0.5 0.6 time (s) d 1 , d 2 Λ = 1, d2 Λ = 2, d2 Λ = 4, d2 Λ = 1, d1 Λ = 2, d1 Λ = 4, d1 (b) Traditional Method. Figure 4.3: Method Comparison for Example 1. 43 Chapter 4. Results 0 1 2 3 4 5 0 0.5 1 1.5 2 2.5 3 3.5 time (s) d 1 , d 2 τ = 0.33, d2 τ = 1, d2 τ = 2, d2 τ = 0.33, d1 τ = 1, d1 τ = 2, d1 (a) New Method. 0 1 2 3 4 5 0 0.5 1 1.5 2 2.5 3 3.5 time (s) d 1 , d 2 Λ = 1, d2 Λ = 8, d2 Λ = 16, d2 Λ = 1, d1 Λ = 8, d2 Λ = 16, d2 (b) Traditional Method. Figure 4.4: Method Comparison for Example 2. 44 Chapter 5 Conclusions We have presented a new method for animating deformable objects. This new method is inspired by physically-based methods, but is adapted to critically damp all modes and to provide an intuitive interface for animators. A new and much simpler form of control over the damping properties of these deformable materials is achieved. The method is capable of producing plausible motion, even when the requested action is not physically possible. The new framework removes the need to choose unintuitive damping parameters, while still allowing for inertial effects and ballistic motion, pro- viding an interesting middle-ground between quasi-static simulation and full physical simulation. The method maintains temporal coherence, and uses a framework compatible with techniques for simulating other materials. 5.1 Future Work There are several advances this method could take. As well, further analysis could lead to a better understanding of traditional damping, and perhaps more intuitive control over traditional damping parameters for a physically accurate simulation. 45 Chapter 5. Conclusions 5.1.1 Oscillation Control In some situations, possibly for artistic effect, the oscillating behaviour of an under-damped system is desired. In order to incorporate this into the new method, we phrase a Backward Euler step of elastic forces in the minimiza- tion framework. Since elastic force is the negative gradient of the elastic potential energy, this step can be written as vnew = v +∆tM−1 (−∇PE(xnew)) (5.1) which solves min vnew 1 2 (∆v)TM(∆v) + PE(xnew) (5.2) where ∆v = vnew − v and xnew = x + ∆tvnew. The new velocity is a function of stress, and therefore so is ∆v. We can then define the new kinetic energy-like term as KEosc(σ) = 1 2 ( ∆v(σ) )TM(∆v(σ)) (5.3) In order to combine the methods, we can blend between the fully elastic term and the kinetic energy in (3.10), min ∫ (κ)KEosc + (1− κ)KE + τ∆tPE (5.4) With κ = 0, we have the new method, (2.22), and with κ = 1, the system is no longer damped. For 0 < κ < 1, we have a method with simple control over oscillations. 46 Chapter 5. Conclusions 5.1.2 Plasticity An important extension for this work is incorporating plastic flow. An object undergoes plastic deformation when the stress, σ, reaches the yield limit, σy. With the current framework, the modified total energy is minimized with respect to unconstrained stress. By adding constraints on stress, the calculated deformation would be purely elastic, and when the yield surface is reached, additional plastic deformation could be calculated. min σ |σ|≤σy ∫ Ω KE(σ) + τ ∆t PE(σ) (5.5) The complications arise when the constraint becomes active and the object permanently changes, affecting even the material coordinates. The minimization may need to be completely restarted when this happens, and combined with the fact that permanent deformations remove some opportu- nities for pre-computations, the method may be quite slow without creative optimization. 5.1.3 Extending the Energy Minimization Framework Since the underlying energy optimization structure can be applied to many other physical phenomena, from frictional contact to compressible flow, this new method lends itself to coupling with other simulations. For example, frictional contact between controlled-damping cloth and flesh could be calcu- lated as a minimization of total kinetic energy and scaled potential energies with respect to internal stresses and contact forces. The contact forces would be subject to inequality constraints, requiring the normal component to be 47 Chapter 5. Conclusions positive and the tangential components to be within the friction cone. This could be further immersed in, and coupled to, an incompressible fluid sim- ulation, by adding the kinetic energy of the fluid to the objective function, and adding fluid pressure as an independent variable. 48 Bibliography [1] Christopher Batty, Florence Bertails, and Robert Bridson. A fast varia- tional framework for accurate solid-fluid coupling. In SIGGRAPH ’07: ACM SIGGRAPH 2007 papers, page 100, New York, NY, USA, 2007. ACM Press. [2] David E. Breen, Donald H. House, and Michael J. Wozny. Predicting the drape of woven cloth using interacting particles. In SIGGRAPH ’94: Proceedings of the 21st annual conference on Computer graphics and interactive techniques, pages 365–372, New York, NY, USA, 1994. ACM Press. [3] R. Bridson, S. Marino, and R. Fedkiw. Simulation of clothing with folds and wrinkles. In SIGGRAPH ’05: ACM SIGGRAPH 2005 Courses, page 3, New York, NY, USA, 2005. ACM Press. [4] P. Dorato, C. Abdallah, and V. Cerone. Linear-Quadratic Control: An Introduction. Prentice-Hall, Englewood Cliffs, NJ, USA, 1995. [5] O. Etzmuss, J. Gross, and W. Strasser. Deriving a particle system from continuum mechanics for the animation of deformable objects. 49 Bibliography Transactions on Visualization and Computer Graphics, 9(4):538–550, 2003. [6] Allen Van Gelder. Approximate simulation of elastic membranes by triangulated spring meshes. J. Graph. Tools, 3(2):21–42, 1998. [7] S. Gibson and B. Mirtich. A survey of deformable modeling in com- puter graphics. Technical Report TR-97-19, Mitsubishi Electric Re- search Lab., 1997. [8] Gentaro Hirota, Susan Fisher, Andrei State, Chris Lee, and Henry Fuchs. An implicit finite element method for elastic solids in contact. In Proc. of Computer Animation, pages 136–146, 2001. [9] Mary K. Kaiser, Dennis R. Proffitt, Susan M. Whelan, and Heiko Hecht. Influence of animation on dynamical judgments. Journal of Experimen- tal Psychology: Human Perception and Performance, 18(3):669–690, 1992. [10] L. Kharevych, Weiwei Yang, Y. Tong, E. Kanso, J. E. Marsden, P. Schröder, and M. Desbrun. Geometric, variational integrators for computer animation. In SCA ’06: Proceedings of the 2006 ACM SIG- GRAPH/Eurographics symposium on Computer animation, pages 43– 51, Aire-la-Ville, Switzerland, Switzerland, 2006. Eurographics Associ- ation. [11] Victor J. Milenkovic and Harald Schmidl. Optimization-based anima- tion. In SIGGRAPH ’01: Proceedings of the 28th annual conference on 50 Bibliography Computer graphics and interactive techniques, pages 37–46, New York, NY, USA, 2001. ACM Press. [12] Gavin S. P. Miller. The motion dynamics of snakes and worms. In SIGGRAPH ’88: Proceedings of the 15th annual conference on Com- puter graphics and interactive techniques, pages 169–173, New York, NY, USA, 1988. ACM Press. [13] Andrew Nealen, Matthias Mller, Richard Keiser, Eddy Boxerman, and Mark Carlson. Physically based deformable models in computer graph- ics. Eurographics State of the Art Report, 2005. [14] Jorge Nocedal and Stephen J. Wright. Numerical Optimization. Springer, 1999. [15] Manfred Nusseck, Julien Lagarde, Benoit Bardy, Roland Fleming, and Heinrich H. Bülthoff. Perception and prediction of simple object inter- actions. In APGV ’07: Proceedings of the 4th symposium on Applied perception in graphics and visualization, pages 27–34, New York, NY, USA, 2007. ACM Press. [16] James F. O’Brien, Adam W. Bargteil, and Jessica K. Hodgins. Graph- ical modeling and animation of ductile fracture. In SIGGRAPH ’02: Proceedings of the 29th annual conference on Computer graphics and interactive techniques, pages 291–294, New York, NY, USA, 2002. ACM Press. [17] James F. O’Brien and Jessica K. Hodgins. Graphical modeling and an- imation of brittle fracture. In SIGGRAPH ’99: Proceedings of the 26th 51 Bibliography annual conference on Computer graphics and interactive techniques, pages 137–146, New York, NY, USA, 1999. ACM Press/Addison-Wesley Publishing Co. [18] Hidehiko Okabe, Haruki Imaoka, Takako Tomiha, and Haruo Niwaya. Three dimensional apparel cad system. In SIGGRAPH ’92: Proceedings of the 19th annual conference on Computer graphics and interactive techniques, pages 105–110, New York, NY, USA, 1992. ACM Press. [19] Carol O’Sullivan, John Dingliana, Thanh Giang, and Mary K. Kaiser. Evaluating the visual fidelity of physically based animations. In SIG- GRAPH ’03: ACM SIGGRAPH 2003 Papers, pages 527–536, New York, NY, USA, 2003. ACM Press. [20] A. Pentland and J. Williams. Good vibrations: model dynamics for graphics and animation. In SIGGRAPH ’89: Proceedings of the 16th annual conference on Computer graphics and interactive techniques, pages 215–222, New York, NY, USA, 1989. ACM. [21] Stephen M. Platt and Norman I. Badler. Animating facial expressions. In SIGGRAPH ’81: Proceedings of the 8th annual conference on Com- puter graphics and interactive techniques, pages 245–252, New York, NY, USA, 1981. ACM Press. [22] Dennis R. Proffitt and David L. Gilden. Understanding natural dy- namics. Journal of Experimental Psychology: Human Perception and Performance, 15(2):384–393, 1989. 52 Bibliography [23] J. Teran, S. Blemker, V. Ng Thow Hing, and R. Fedkiw. Finite volume methods for the simulation of skeletal muscle. In SCA ’03: Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 68–74, Aire-la-Ville, Switzerland, Switzerland, 2003. Eurographics Association. [24] Joseph Teran, Eftychios Sifakis, Geoffrey Irving, and Ronald Fedkiw. Robust quasistatic finite elements and flesh simulation. In SCA ’05: Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 181–190, New York, NY, USA, 2005. ACM Press. [25] Xiaoyuan Tu and Demetri Terzopoulos. Artificial fishes: physics, lo- comotion, perception, behavior. In SIGGRAPH ’94: Proceedings of the 21st annual conference on Computer graphics and interactive tech- niques, pages 43–50, New York, NY, USA, 1994. ACM Press. [26] Keith Waters. A muscle model for animation three-dimensional facial expression. In SIGGRAPH ’87: Proceedings of the 14th annual con- ference on Computer graphics and interactive techniques, pages 17–24, New York, NY, USA, 1987. ACM Press. 53
The Open Collections website will be unavailable July 27 from 2100-2200 PST ahead of planned usability and performance enhancements on July 28. More information here.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Controllable, non-oscillatory damping for deformable...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Controllable, non-oscillatory damping for deformable objects Young, Herbert David 2007
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Controllable, non-oscillatory damping for deformable objects |
Creator |
Young, Herbert David |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | This thesis presents a new method for the controllable damping of deformable objects. The method evolves from physically based techniques; however, it allows for non-physical, but visually plausible motion. This flexibility leads to a simple interface, with intuitive control over the behaviour of the material. This method is particularly suited for strongly damped materials, which account for the majority of objects of interest to animation, since it produces non-oscillatory behaviour. This is similar to critical damping, except that it affects all modes independently. The new method is based on the minimization of a slightly modified version of total energy. This framework can be used to simulate many other physical phenomena, and therefore lends itself to coupling with other simulations. Implementation details for a simple example are given. Results are shown for varying parameters and compared to those produced by a traditional method. |
Extent | 392369 bytes |
Subject |
damping elastic deformable energy minimization |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2007-12-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0051414 |
URI | http://hdl.handle.net/2429/217 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2008-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2008_spring_young_herbert.pdf [ 383.17kB ]
- Metadata
- JSON: 24-1.0051414.json
- JSON-LD: 24-1.0051414-ld.json
- RDF/XML (Pretty): 24-1.0051414-rdf.xml
- RDF/JSON: 24-1.0051414-rdf.json
- Turtle: 24-1.0051414-turtle.txt
- N-Triples: 24-1.0051414-rdf-ntriples.txt
- Original Record: 24-1.0051414-source.json
- Full Text
- 24-1.0051414-fulltext.txt
- Citation
- 24-1.0051414.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0051414/manifest