Eulerian Finite Volume Method for Musculoskeletal Simulation and Data-driven Activation by Ye Fan B.Eng. Shanghai Jiaotong University, 2010 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Computer Science) The University Of British Columbia (Vancouver) June 2013 c Ye Fan, 2013 Abstract This thesis describes a solid simulation method and its application to musculoskeletal simulation. The presented solid simulation method features Eulerian discretization and avoids mesh tangling during large deformation. Unlike existing Eulerian solid simulation methods, our method applies to elastoplastic material and volume-preserving material. To further increase the utility of Eulerian simulations for solids, we introduce Lagrangian modes to the simulation and present a new solver that handles close contact while simultaneously distributing motion between the Lagrangian and Eulerian modes. This Eulerian-on-Lagrangian method enables unbounded simulation domains and reduces the time step restrictions that often plague Eulerian simulation. We also introduce a framework for simulating the dynamics of musculoskeletal systems, with volumetric muscles and a novel muscle activation model. Muscles are simulated using the solid simulator developed and therefore enjoys volume preservation which is crucial for accurately capturing the dynamics of muscles and other biological tissues. Unlike previous work, in our system muscle deformation is tightly coupled to the dynamics of the skeletal system, and not added as an after effect. Our physiologically based muscle activation model utilizes knowledge of the active shapes of muscles, which can be manually drawn or easily obtained from medical imaging data. Finally we demonstrate results with models derived from MRI data and models designed for artistic effect. ii Preface All the work presented in this thesis were guided by Dr. Dinesh K. Pai and were described in co-authored research papers: • Ye Fan, Joshua Litven, David I.W. Levin, and Dinesh K. Pai, Eulerian-onLagrangian Simulation (ACM Transactions on Graphics, June 2013) • Ye Fan, Joshua Litven, and Dinesh K. Pai, Volumetric Musculoskeletal Simulation and Data-driven Activation, 2013 (in preparation) The work presented in Chapter 4 has resulted in the first paper and will be presented at SIGGRAPH 2013. Dr. David I.W. Levin originally designed and developed the Eulerian Solid simulation method. After that, we have greatly extended his work in various ways and made the method more practical. Dr. Levin provided theoretical guidance and I implemented the Eulerian-on-Lagrangian simulation algorithm with GPU acceleration. Joshua Litven designed and implemented the Quadratic Programming solver. Chapter 5 and Section 3.3 were mainly based on the second paper, where I further extended Eulerian-on-Lagrangian simulation’s applicability to volume preserving solids and applied it to musculoskeletal simulation. Additionally, Dr. Pai and I designed and implemented a novel muscle activation model. In this joint work, Joshua Litven contributed his expertise of developing an iterative QP solver used for resolving incompressibility and contact. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . 3 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 Simulating Deformable Objects . . . . . . . . . . . . . . . . . . 4 2.2 Modeling Musculoskeletons . . . . . . . . . . . . . . . . . . . . 6 Finite Volume Method for Eulerian Solid Simulation . . . . . . . . . 8 3.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Volume Preservation . . . . . . . . . . . . . . . . . . . . . . . . 13 3.4 Contact and Collision . . . . . . . . . . . . . . . . . . . . . . . . 20 3.4.1 21 2 3 Surface Reconstruction . . . . . . . . . . . . . . . . . . . iv 3.4.2 Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Eulerian-on-Lagrangian Simulation . . . . . . . . . . . . . . . . . . 25 4.1 Rigid Motion Eulerian-on-Lagrangian Simulation . . . . . . . . . 27 4.2 Linear Modal Eulerian-on-Lagrangian Simulation . . . . . . . . . 29 4.3 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . 29 4.4 Exploiting Redundancies . . . . . . . . . . . . . . . . . . . . . . 32 4.5 Solving the Dynamics in the Presence of Contact . . . . . . . . . 35 Musculoskeletal Simulation and Activation . . . . . . . . . . . . . . 38 5.1 Bones and Tendons . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.2 Coupling of Bones, Tendons and Muscles . . . . . . . . . . . . . 42 5.3 Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.4 Data-driven Muscle Activation . . . . . . . . . . . . . . . . . . . 46 Results and Examples . . . . . . . . . . . . . . . . . . . . . . . . . . 49 6.1 Volume-Preservation . . . . . . . . . . . . . . . . . . . . . . . . 49 6.2 Eulerian-on-Lagrangian Simulation . . . . . . . . . . . . . . . . 52 6.3 Musculoskeletons . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.3.1 Eye muscles . . . . . . . . . . . . . . . . . . . . . . . . 55 6.3.2 Arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . 60 7.1 Limitations and Future Work . . . . . . . . . . . . . . . . . . . . 60 7.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.5 4 5 6 7 v List of Tables Table 4.1 Notation used . . . . . . . . . . . . . . . . . . . . . . . . . . Table 6.1 For each example, the size of the spatial grid (Dim), maximum 26 GPU memory used, total number of time steps taken, and average runtimes in milliseconds for each stage of the algorithm is shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 6.2 54 Statistics for a typical time step of the iterative QP solver for varying dimensions of the human arm example. For each resolution, The number of DOFs and constraints of each type is shown, as well as runtimes of applying the KKT linear operator and preconditioner, in milliseconds. . . . . . . . . . . . . . . . vi 54 List of Figures Figure 3.1 Rectangular simulation grid . . . . . . . . . . . . . . . . . . Figure 3.2 a) An object with volume Ω inside of the collocated grid b) 12 Each grid cell is decomposed into subcells Ωi j to transform flux surface integrals along Γc to integrals along the inner subcell faces Γa and Γb , depicted as bold lines. . . . . . . . . . . 17 Figure 3.3 Two plastic rods are colliding with a rigid rod. . . . . . . . . 24 Figure 4.1 The continuous domains used for the Eulerian-on-Lagrangian simulator as well as the discrete grid structure stores in the Eulerian domain. The mapping between the domains is described using a set of generalized configuration variables. . . . . . . . Figure 4.2 26 The block structure of the Eulerian-on-Lagrangian mass matrix using rigid motion. The diagonal blocks are Lagrangian and Eulerian respectively while the off-diagonal blocks represent the coupling terms of the system. . . . . . . . . . . . . . Figure 4.3 30 Due to redundancy of the Lagrangian modes, the motion of an object can be specified by advection through the fixed Eulerian grid (top right), or fixing the object in the grid and moving the grid itself via the Lagrangian degrees of freedom (bottom right). 34 Figure 5.1 Muscle, tendon and bone structure . . . . . . . . . . . . . . . Figure 5.2 A 2D rigid body system with an inextensible tendon (orange line) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii 38 41 Figure 5.3 The mapping between active space, material space and physical space . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Figure 5.4 Muscle deformation by manually specified activation shape . . 47 Figure 6.1 The figure shows the volume change for three simulations in comparison. y-axis: object volume in percentage; x-axis: total simulation time in seconds. The red line is for the non-volumepreserving solid while the green and blue lines are for volume preserving solids. The green one is based on ∇·v = 0 condition and the blue one is based on our method . . . . . . . . . . . . Figure 6.2 50 The figure shows the volume change for objects with different stiffness. The blue line representing our method is for reference. The object’s stifness is 1000kP. The other three lines represent objects with Young’s Modulus 1000kP, 1250kP and 1500kP respectively, simulated without stabilization . . . . . Figure 6.3 51 Comparison of volume-preserving and non-volume-preserving bunnies dropping to the ground. Left: volume-preserving bunny. Right: purely elastic bunny . . . . . . . . . . . . . . . . . . . Figure 6.4 Left: Purely Eulerian motion. Middle: Purely Lagrangian motion. Right: Overall motion . . . . . . . . . . . . . . . . . . . Figure 6.5 55 Top Row: Purely Eulerian translation. Bottom Row: Eulerianon-Lagrangian Translation . . . . . . . . . . . . . . . . . . . Figure 6.6 52 55 The angular velocity of an Eulerian-on-Lagrangian object (Eon-L) and a purely Eulerian object undergoing rotation. Both were given an initial angular velocity of π radians/s. . . . . . Figure 6.7 56 The time step, in seconds, of an Eulerian-on-Lagrangian simulation. Note that during contact the time step is automatically reduced. After contact the time step increases until it reaches the rendering frame rate. . . . . . . . . . . . . . . . . . . . . Figure 6.8 56 The timestep size for both the purely Eulerian and Eulerian-onLagrangian simulations for the two ball collision. Note that in most cases the EoL time step is an order of magnitude greater than that of the purely Eulerian simulation. . . . . . . . . . . viii 57 Figure 6.9 ˜ middle:the Left: the MRI data of the active configuration X; reference configuration X; right: simulations of the activated muscles in physical space x . . . . . . . . . . . . . . . . . . . 58 Figure 6.10 Up: the upper arm in the rest shape; Down: the Biceps and the Brachialis are activated . . . . . . . . . . . . . . . . . . . . . ix 59 Acknowledgments First, I’m giving my deep gratitude to my supervisor Dr. Dinesh K. Pai, who guided me to mechanics and offered numerous opportunities for letting me do my dream work. Also I would like to thank Dr. Robert Bridson for reading my thesis and providing invaluable advices. I would like to extend my acknowledgements to all the colleagues and friends at the Sensorimotor Lab at UBC, especially to Dr. David I.W. Levin for the enormous help in research and life and for teaching me everything about Eulerian solids. Special thanks to Joshua Alexander Litven, who is my great friend, collaborator and English teacher. I also appreciate all the help from Dr. Bin Wang, who coded collision detection for the Eulerian-on-Lagrangian simulation. Finally, I sincerely thank my family for all the love and support. x Chapter 1 Introduction This thesis concerns simulation of elastic objects with the ultimate goal of musculoskeletal simulation. Though abundant related research results exist in different areas, we aim at graphics applications where efficiency, robustness and visual effects are given higher priority to the accuracy of the methods. In graphics, simulating deforming objects in a physically based way has drawn much attention. Recently, Levin et al. [2011b] proposed an Eulerian solid method that has proven to be effective for simulating soft objects in large deformation and close contact in an easy and robust way. This makes the Eulerian method an ideal candidate for muscle simulation because muscles undergo large deformation and are always in contact with neighboring tissues, such as other muscles, fat and skin. However, the Eulerian method suffers from many shortcomings that are in fact its defining characteristics, for example, the fixed simulation domain and small time steps. In this thesis, we have improved this method in two major aspects and significantly increased the utility of Eulerian simulations for materials other than fluids. The first contribution is Eulerian-on-Lagrangian simulation method, which alleviates many of the shortcomings of purely Eulerian and purely Lagrangian methods. It features unbounded simulation domains, larger time steps and less numerical dissipation. The other contribution is extending the applicability of the method to elasto-plastic material and incompressible material. For elasto-plastic material, the Eulerian method avoids artificial plasticity because it doesn’t need constant remeshing. Due to such an advantage, we have further proposed a muscle acti1 vation method that follows the same approach of simulating elasto-plasticity. The ability to simulate volume-preserving solids is equally important since it allows the method to simulate biological tissues which are mostly composed of water, such as muscles. The second contribution serves as the building block for the second part of the thesis: musculoskeletal simulation and activation. From the perspective of graphics research, understanding the muscle mechanisms and reproducing interesting muscle behavior have obsessed generations of talented scientists. To highlight a few fascinating results, Sueda et al. [2008] has successfully modeled human hand and forearm, which is one of the most expressive parts of human body. Lee et al. [2009] modeled the full human upper body using line muscles for activation and volumetric muscles for visualization. In the second part of the thesis, we focus on the problem of muscle deformation and activation, which is indispensible for realistic and artistic character animation. Imagine Popeye the sailor, who bulges his biceps in an exaggerated manner upon consuming a healthy dose of spinach. If we want to redo such a scenario in a biologically reasonable way, how can we automate the whole process? The presented thesis has a collection of techniques that address such demands. Based on the improved Eulerian solid simulation method described in the first part of the thesis, we have developed a musculoskeletal simulator that handles bones, muscles and tendons in a unified framework. It captures crucial characteristics of the muscles, such as volume preservation and close contact. While many other works exist on passive muscle simulation, they do not have a practical activation model which is crucial for realistic character animation. Admitting that the actual mechanism of muscle activation still remains unknown and an active research topic, we provide a simple but effective method that addresses the needs of graphics. The presented activation model allows modelers to draw artistic shapes for activated muscles and leverages existing biomedical data such as ultra-sound and Magnetic Resonance Imaging (MRI). To be more specific, the novel activation method follows the multiplicative decomposition law that was originally used for simulating elasto-plastic material. As a demonstration, we have utilized actual MRI data for simulating the movement of human eyeballs. 2 1.1 Contributions In summary, the contributions in the thesis include the following • The Eulerian-on-Lagrangian simulation method that alleviates many of the shortcomings of purely Eulerian and purely Lagrangian discretizations. • Extending the Eulerian solid simulation method to incompressible material and elasto-plastic material. • Implementing a musculoskeletal simulator that couples bones, muscles and tendons. • A novel muscle activation model. 1.2 Thesis Organization In Chapter 2, we review the most closely related work in solid simulation and musculoskeletal simulation. Chapter 3 introduces the Eulerian simulation methodology and describes volume preservation and elasto-plasticity in detail. We then discuss Eulerian-on-Lagrangian simulation in Chapter 4. The musculoskeletal simulation framework is described in Chapter 5. Chapter 6 gives the simulation results and Chapter 7 concludes this thesis. 3 Chapter 2 Related Work Simulating the human musculoskeletal system is a fascinating research topic that has attracted efforts from various areas. Our work draws from two distinct research fields: numerics of solid simulation and musculoskeletal modeling. While our motivation is primarily based on a graphical application, these fields are inherently related to this thesis and possess a plethora of results and computational techniques. 2.1 Simulating Deformable Objects Numerical simulation methods for solids have a long and rich history. Especially after the Finite Element Method (FEM) was invented, hyperelasticity simulation had been widely applied in industry for rubber simulation. Since FEM is theoretically well grounded, it is considered as a principled method for numerical simulation. Recently, FEM has equipped itself with modern techniques to address more complicated PDEs. Eulerian methods have become commonplace in computer graphics for fluid simulation. Unlike Lagrangian-type discretization, methods in the Eulerian family take a very different approach. Instead of discretizing the material space, it discretizes the world space. Objects “flow” around in the simulation grid, which makes a lot of sense for fluids. Eulerian methods have also been used to simulate melting solids [Carlson et al., 2002], viscoelastic fluids [Goktekin et al., 2004, Losasso et al., 2006] and elastic or elasto-plastic solids [Banks et al., 2007, Bar- 4 ton and Drikakis, 2010, Kamrin and Nave, 2009, Kamrin et al., 2012, Levin et al., 2011b, Miller and Colella, 2001, Sulsky et al., 1994, Tran and Udaykumar, 2004, Trangenstein, 1994]. For hyperelasticity, Eulerian solid simulation method is pioneered by Kamrin et al. [2009, 2012] and Levin et al. [2011b]. Their work has shown great potential for simulating soft objects under large deformation using a reference map technique. Based on their approach, we develop a general simulator for incompressible solids, which extends the method’s modeling capabilities to that of general materials. Typically, resistance to volume change is approximated by tuning the Poisson ratio of a material. This is normally effective for virtual objects. However, for certain materials such as biological tissue, accurate volume preservation is crucial. As the Poisson ratio approaches a value of ≈ 0.5, it induces a numerical difficulty known as the “locking phenomenon”. The classical solution is through mixed FEM formulation [Hughes, 2000] with carefully chosen elements, or a stabilization term [Hughes, 2000, Misztal et al., 2012] is introduced. Recent work by Irving et al. [2007] and Patterson et al.[2012] belong to these two approaches as well. Unlike the quasi-incompressibility method[Patterson et al., 2012], our work linearizes the incompressibility condition on the velocity level and features dynamic simulation. As for stabilization, we take the same approach as the method presented by Misztal et al. [2012] However, Eulerian simulations suffer from several inherent drawbacks such as spatio-temporal time step limitations [Osher and Fedkiw, 2002], dissipation caused by certain advection schemes [Stam, 1999] and the necessity of defining a computational domain that extends throughout the simulation domain [Foster and Metaxas, 1996]. For Eulerian fluid simulations these shortcomings have been addressed by using more complicated, conservative advection schemes [Lentine et al., 2011] and data structures [Wang et al., 2005]. An alternate approach involves replacing the Eulerian computational machinery with Lagrangian analogs [Premoˇze et al., 2003]. However in these cases we are faced with losing the advantages conferred by Eulerian simulation. Methods which maintain Eulerian and Lagrangian characteristics simultaneously can prove to be useful alternatives to these either-or approaches. Such methods have been used effectively for simulation of constrained, one-dimensional strands [Sueda et al., 2011], dissipation free advection [Zhu and Bridson, 2005], robust pressure-projection in particle-based 5 fluid simulations [Narain et al., 2010, Raveendran et al., 2011] and as a means of solid fluid coupling in finite-element simulations [Belytschko and Kennedy, 1978]. Additionally, methods which translate [Shah et al., 2004] or rotate and scale [Poludnenko and Khokhlov, 2007] the simulation domain have been used in fluid simulation and flexible objects [Terzopoulos and Witkin, 1988]. We will take a different approach from the methods listed above. Instead of assigning algorithmic steps or discrete points to be Lagrangian or Eulerian we will imbue generalized configuration variables with these characteristics via a sequence of mappings. This will allow us to assign Lagrangian and Eulerian behavior to the different modes of a motion (not just to the nodes of a domain). We will give two specific examples of Eulerian-on-Lagrangian methods but note that the algorithm described allows a completely general decomposition to be applied. As opposed to the Arbitrary-Lagrangian-Eulerian method [Belytschko et al., 2000] (ALE) and Lagrangian methods which use remeshing to avoid ill-conditioned elements [Bargteil et al., 2007, Wicke et al., 2010], our method describes a principled numerical “glue” that can be used to merge Eulerian and Lagrangian Simulations, both described using arbitrary generalized coordinates, and inherit the advantages of both schemes. 2.2 Modeling Musculoskeletons Simulating musculoskeletal systems, with 3D volumetric muscles, has been an active research area in a multitude of disciplines, from biomechanics [Blemker and Delp, 2005] to computer animation [Lee et al., 2012]. Within graphics, Chen et al. [1992] first introduced muscle simulation to character animation. To highlight a few recent results, Teran et al. [2003, 2005] build and simulate skeletal muscles from the Visible Human Data Set; Lee et al. [2009] simulate the entire human upper body using FEM and the line muscle model; Sueda et al. [2008] focus on the dynamics of thin musculotendons and simulate an anatomy-based hand. For modeling volumetric muscles on the macroscopic level, the fiber reinforced model is commonly adopted in graphics and biomechanics [Teran et al., 2003] using “Hill-type” models of muscle activation. Zajac et al. [1989] provide an extensive review of muscle and tendon models that are commonly used in biomechanics. 6 However, the actual properties of muscles [Nishikawa et al., 2012] and their fiber architecture [Levin et al., 2011a] are not sufficiently well known and are difficult to measure. Alternatively we focus on developing an activation model based on what we can obtain more easily, either from medical images or from an artist. This allows us to use simpler hyperelastic constitutive laws since we do not rely on fiber information to activate the muscles. Our activation model works by introducing an active state configuration for the material. It shares the same spirit as the method proposed by Nardinocchi et al. [2007] in the sense that both methods follow a multiplicative decomposition of the deformation gradient into a passive part and an active part. The idea originally comes from mechanics of elastoplasticity [Lee, 1968]. Related computational methods appear in graphics for simulating plastic flows [Bargteil et al., 2007, Wicke et al., 2010]. Research in medical image analysis is a good complement to our muscle activation model since there exist advanced methods for obtaining shapes of active muscles from ultrasound and MRI [Levin et al., 2011a, Wei and Pai, 2008]. Those data can be readily used for active state configurations. 7 Chapter 3 Finite Volume Method for Eulerian Solid Simulation We have briefly mentioned the advantages of simulating muscles using the Eulerian solids method. In this section, we give the background of the method and discuss volume preservation in detail. First presented by Levin et al. [2011b], the Eulerian solids method addresses the fundamental difficulties of simulating solids under large deformation and close contact. In previous work on solid simulation, the prevailing methods are based on Lagrangian discretization, which is ideal for solids such as metals. As for softer materials, biological tissues in particular, those Lagrangian simulation methods easily suffer from mesh tangling during large deformation. If mesh quality becomes too bad to produce meaningful results, remeshing is required. Unfortunately remeshing algorithms are nontrivial to implement. On the other hand, Eulerian methods discretize the physical space, providing a different solid simulation paradigm that does not require sophisticated remeshing. Due to this nice property, different objects in contact can share one common grid cell for contact formulations. The corresponding Lagrangian methods need sophisticated ways to constrain the model, for example, the mortar method [Belgacem et al., 1998, Hild, 2000]. In this sense, Eulerian methods are advantageous for simulating multiple muscles in consistent contact. In this thesis, we describe a muscle simulator based on the Eulerian framework. To accommodate particular requirements, such as vol8 ume preservation and skeletal couplings, we greatly extend the original work in various ways. The new contributions are highlighted in later sections. First of all, to make the thesis self-contained, the fundamental concepts and technical details are given below. 3.1 Kinematics Kinematics is the study of motion and deformation without reference to the cause. For the Eulerian solid formulation, we are working with the material space X and physical space x. The motion of an object is defined as a mapping between X and x. As long as objects do not interpenetrate each other, the mapping is continuous and we define it as x = φ (X), which implies a one-to-one correspondence between the points in these two spaces. Inversely, X = φ −1 (x) uniquely exists. For spatial discretization, we either store the forward mapping ψ or the inverse mapping φ −1 . The Eulerian method takes the latter approach. It discretizes x on a fixed grid and stores the corresponding positions in X on each of the grid nodes. This way, the Eulerian method defines the discretized inverse mapping φ¯ −1 . The discretized forward mapping φ¯ implicitly exists. For notational simplicity, we follow the common practice and use X to denote the mapping φ −1 . The inverse deformation gradient F−1 is defined as F−1 = ∂X ∂x (3.1) Computing F is important because it measures object deformation and is indispensible for evaluating strain and stress. In the Eulerian setting, F is easily computed by inverting the inverse deformation gradient F= ∂x = ∂X ∂X ∂x −1 (3.2) The determinant of the deformation gradient, J = det(F), measures the volume change of an infinitesimal volume element in material space. It is given as dv = JdV (3.3) where dv is the volume in physical space and dV is the volume in material space. 9 3.2 Dynamics We formulate object motion by following Gauss’s principle of least constraints [Lanczos, 1986]. It minimizes the Gibbs’ function while satisfying all the constraints. The constrained optimization problem is given as 1 T a Ma − fT a 2 subject to a ∈ A , minimize (3.4) where M is a mass matrix, a and f are the accelerations and forces defined in physical space, and A is the set of accelerations that are consistant with the constraints; these constraints ensure volume preservation and non-interpenetration between objects. f includes body forces, elastic forces and inertial coupling forces. Elastic forces are computed through a finite volume approach. More details are available in the paper by Levin et al. [2011b]. Contact forces and hydrostatic forces are implicitly calculated by solving the resulting KKT system. In our formulation, velocities are selected as the primary Degrees of Freedoms (DOFs). We discretize acceleration a in time For the advective term an+1 = vn+1 − vn ∂ vn n+1 + v ∆t ∂x ∂ vn n+1 , ∂x v we use the velocity vn at the current time step (3.5) in order to make the cost function in Equation 3.4 quadratic. Substituting Equation 3.5 into Equation 3.4, we get vn+1 by solving the quadratic programming (QP) problem vn+1 =argmin v 1 T v (I + ∆t∇vn )T M(I + ∆t∇vn )v − (∆tf + Mvn )T v 2 (3.6) subject to GC v ≥ 0 GV v = c, where GC and GV represent the contact and the volume-preserving constraint matrices repectively. c represents possible non-zero values of the volume-preserving constraint on the right hand side. ∂v ∂xv can also be explicitly discretized using the 10 current vn as ∂v ∂ vn n v= v ∂x ∂x (3.7) Substituing Equation 3.7 into the definition of a, Equation 3.5 becomes an+1 = vn+1 − vn ∂ vn n + v ∆t ∂x (3.8) Substituting Equation 3.8 into Equation 3.4 gives an alternative optimization problem for solving for vn+1 1 T v Mv − (∆tf + M(vn − ∆t(vn · ∇)vn ))T v 2 vn+1 =argmin v (3.9) subject to GC v ≥ 0 GV v = c, The advective term (v · ∇)v does not appear in the time discretization of Lagrangian methods. Due to the inherent non-linearity, solving it both accurately and efficiently is nontrivial. In graphics, the practical approach is time splitting. The basic idea is to split a complicated equation into several parts and solve them separately. More explanations and examples about time splitting are well described in the book by Bridson [2008]. We formulate the dynamics based on Equation 3.9, but the implementation is similar to splitting. First, we evaluate the advective term using the semi-Lagrangian method, which is unconditionally stable. v∗ = vn − ∆t(vn · ∇)vn (3.10) Then, we get vn+1 by solving vn+1 =argmin v 1 T v Mv − (∆tf + Mv∗ )T v 2 subject to GC v ≥ 0 GV v = c, where v∗ is already known from Equation 3.10. 11 (3.11) After vn+1 is obtained, X is updated by solving ˙ = ∂X + ∂Xv = 0 X ∂t ∂x (3.12) which bears the physical meaning that the object in material space is always fixed. Equation 3.12 can also be thought of as advecting material coordinates X in physical space. Algorithm 1 A high level overview of the dynamics in one time step. 1: 2: 3: 4: 5: v∗ = advect(vn ) Compute all the force terms f Formulate the constraints Solve Equation 3.11 for vn+1 Solve Equation 3.12 for X The steps of updating v and X in one time step are given in the Algorithm 1 So far, we have explained the dynamics and the time discretization in the algo- rithm, which plays a crucial role for the robustness and the efficiency of the simulation. Implementing a more sophisticated time integration scheme is possible, but involves large amount of efforts. We will leave it as future work. Figure 3.1: Rectangular simulation grid 12 Next, we describe the spatial discretization for the Eulerian Finite Volume Method. Figure 3.1 illustrates the spatial discretization in our implementation where a structured rectangular grid is used. We chose a rectangular grid due to ease of implementation and memory efficiency, but the Eulerian method itself does not exclude other popular discretizations like tetrahedral meshes. For the discretization scheme, variables X are stored on the grid nodes (solid circles in Figure 3.1). They define a mapping from physical space to material space. Velocities v are stored at the same places as X. Deformation gradient F, strain and stress are computed and stored at the nodes (solid triangles) of the dual grid (the red dotted grid in Figure 3.1). The mass matrix M is computed in a standard manner. Like applications in graphics and haptics that focus on good performance rather than accuracy, we lump the mass matrix to the diagonal. Such diagonal structure has been well utilized for the solver implementation. 3.3 Volume Preservation Muscles, being biological tissue, are mostly water and therefore nearly incompressible. Capturing volume preservation in muscles is important for interesting behavior such as flexing and bulging. Ng-Thow-Hing et al. [2001] simulates muscles by approximating this constraint using spring-like forces based on volume change. Others follow the quasi-incompressibility method, originally proposed by Simo et al. [1991], with finite element methods for muscle simulation [Teran et al., 2005, Weiss et al., 1996]. The problem of this formulation is that volume-preservation is not guaranteed after volume-preserving forces are applied, especially in the presence of other constraints, such as contact. Stiffness of such forces needs manual tuning and the process can be tedious. Alternatively, since material space stores the initial configuration of the object, the volume-preserving constraint can be formulated by the hard constraint J = 1. Unfortunately, in our simulation framework described in Section 3.2, J = 1 is inconvenient to work with because a) v is selected as the unknowns and all other constraints are formulated on the velocity level. b) J = 1 is non-linear, but our QP solver is only optimized for linear constraints. Therefore, we linearize J = 1 by taking the time derivative and take the fluid ap- 13 proach to solid simulation similar to the method by Irving et al. [2007]. The method is compatible with any constituitive model. In fluid simulation, incompressibility constraint is formulated as ∇ · v = 0 and solved by pressure projection. The pressure projection step decomposes the velocity field into a divergence-free part and a curl-free part (also known as the Hodge decomposition) by solving a Poisson equation for the hydrostatic pressures. Notice that ∇ · v = 0 applies on the velocity level. Local drift of the volume can easily accumulate. The method by Irving et al.[2007] avoids these errors by post-stabilization, which projects the positions to exactly conserve volume. In our work, we present a new stabilization technique that is simple to implement and exhibits good long term behavior. Typically in fluid, this problem is avoided since there is no position variable. But fluid simulations are not totally immune from this issue. If level set methods are used for capturing the surface, liquids may smear out over long simulations. This is due to a known issue of level set method: mass is not conserved neither locally nor globally. To combat such volume loss, effective methods [Chentanez and M´uller, 2012] have been proposed in graphics. For our Eulerian solid simulation, we do not need to worry about these errors caused by implicit surface capturing because of the reference mapping technique. In this section, we briefly give the formulation for volume-preserving solids based on ∇ · v = 0 and details of the stabilization technique. Following the notations in Section 3.1, the motion of an object is defined as the mapping φ −1 (x) = X. (3.13) Recalling from Section 3.1, the determinant J of the deformation gradient F defines the volume change as dv = JdV (3.14) Volume-preservation indicates that for a chosen piece of material, the volume in the material space is equal to the volume in the physical space, dV = dv. Therefore, J=1 (3.15) Equation 3.15 looks rather simple, but is actually non-linear. For computational 14 convenience, in linear elasticity, Equation 3.15 is linearized as J˜ ≈ J + J∇ · u (3.16) where J˜ is the Jacobian after a small change u is made to x. Since in the linear elasticity theory, only very small change is allowed, the incompressibility condition reduces to J∇ · u = 0 (3.17) ∇·u = 0 (3.18) and further The full derivation of Equation 3.16 can be found in the book by Bonet et al. [1997](section 2.3). Here, we only outline a few important steps. From the definition of J and J,˜ the linearization is ∂ (x + u) J˜ = det ∂X ≈ det(F) + Ddet(F)[u] (3.19) where u is the infinitesimal displacement. To proceed, the directional derivative Ddet(F)[u] needs evaluation and has the following expression(see Equation(2.118) in the book by Bonet et al. [1997] for the details) ∂u ∂X (3.20) = Jtr F−1 ∇uF (3.21) = Jtr (∇u) = J∇ · u (3.22) Ddet(F)[u] = Jtr F−1 This linear condition Equation 3.18 is a legitimate simplification for small deformation near the equilibrium state, but it no longer holds physical meaning when the object is under large deformation. Such linearization greatly limits the applicability of the simulator. Instead, we take a different approach and formulate a velocity level constraint by taking the total time derivative on both sides of Equation 3.15 J˙ = 0 15 (3.23) Based on the relationship between time and directional derivatives stated in the book [Bonet and Wood, 1997](Section 3.11.3), the material time derivative of the Jacobi Equation 3.23 can be evaluated as J˙ = DJ[v] (3.24) Recalling the linearization of J in Equation 3.20, Equation 3.24 has the similar derivation and the result is J˙ = J∇ · v (3.25) where J = 0, otherwise, the element is flattened. Therefore, ensuring J˙ = 0 is equivalent to ∇·v = 0 (3.26) Like any other velocity level constraint, Equation 3.26 loses positional information. This is especially an issue for volume-preserving solids because the divergence free velocity field only enforces zero volume change rate. Without compensating the actual volume loss over time, the accumulated numerical drifts can cause permanent volume change to the solid. The overall volume now highly depends on such drifts and is not well predictable. The problem to be solved really comes down to the question “How should we encode positional level information to Equation 3.26 so that numerical drifts can be bounded?”, or “Are there any stabilization methods that can be appropriately borrowed here?” The approach presented here is to enforce the time integral of J˙ be zero, which has the form: tn+1 tn J˙ = 0 J˙+ 0 tn+1 J˙ = 0 (3.27) tn where tn+1 and tn denote the next time step and the current time step respectively. At a typical time step, say tn , it is straightforward to compute J based on X. By substituting Equation 3.25 to Equation 3.27 and assuming J0 = 1, we have Jn − 1 + tn+1 tn 16 J∇ · v = 0 (3.28) To give the expression for the constraint based on Equation 3.28, we need to numerically integrate J∇ · v over time tn and tn+1 . To make the presentation simple, we linearize every term and write the integral as tn+1 J∇ · v = ∆tJn ∇ · vn+1 (3.29) tn Substitute Equation 3.29 back to Equation 3.28, the velocity constraint now is ∇ · vn+1 = 1 − Jn ∆tJn (3.30) Of course, any numerical integration scheme can be employed here, for example tn+1 tn J∇ · v = ∆t (Jn ∇ · vn + Jn+1 ∇ · vn+1 ) 2 (3.31) Jn+1 is a function of Xn+1 which is updated by vn+1 through material derivative, thus Jn+1 is nonlinear to v and its expression is much involved in an Eulerian setting. In practice, we can approximate Jn+1 by 1 and the constraint remains linear on v. Figure 3.2: a) An object with volume Ω inside of the collocated grid b) Each grid cell is decomposed into subcells Ωi j to transform flux surface integrals along Γc to integrals along the inner subcell faces Γa and Γb , depicted as bold lines. 17 From this point on, many classical techniques in Computational Fluid Dynamics (CFD) can be utilized for discretization because the grid is rectangular and fixed. As such, transforming a fluid simulator, e.g. stable fluid [Stam, 1999], into a volume-preserving solid simulator requires few modifications. We now present the details of discretization on a collocated grid using FVM in 2D for illustration purposes; extending to 3D is straightforward. Consider a solid with volume Ω in Figure 3.2. The velocities are stored at the nodes (denoted by solid circles). We are going to formulate the incompressible constraint for each cell. To start, we write the integral of Equation 3.30 over the volume Ω as 1 − Jn dΩ, Ω ∆tJn ∇ · vdΩ = Ω (3.32) Applying the divergence theorem to the left hand side of Equation 3.32 1 − Jn dΩ, Ω ∆tJn ˆ = v · ndΓ Γ (3.33) where Γ = ∂ Ω is the boundary of the volume. nˆ is the outwards pointing normal at the boundary. To evaluate Equation 3.33, we utilize a flux function T, which can be used to evaluate any vector quantity defined on the grid. The numerical evaluation of T is the key to FVM. Generally, the expression of T is TΓ (q) ≡ ˆ q · ndΓ. (3.34) Γ where the subscript Γ denotes the boundary. q is a general vector defined on the grid. Use the flux definition Equation 3.34 and substitute it to Equation 3.33, Equation 3.33 becomes more concise TΓ (v) = 1 − Jn dΩ, Ω ∆tJn (3.35) As for the right hand side of Equation 3.35, we compute the value at the center of Ω and multiply it with the volume, denoted as V . Therefore, Equation 3.35 becomes TΓ (v) = 1 − Jn V ∆tJn 18 (3.36) Now, we only need to evaluate the flux TΓ (v). In order to proceed, we define subcells Ωi j for i, j ∈ {1, 2} as depicted in Figure 3.2, and the flux of the subcells are Ti j (q) = T(q) ∩ Ωi j . (3.37) The total flux is the sum of the fluxes from all the subcells. TΓ (q) = ∑ Ti j (q) (3.38) i, j To compute the flux Ti j (q) we follow the method by Teran et al. [2005]. Using Ωi j ∇ · q constant shape functions implies = 0, so that fluxes may be computed along the inner subcell boundaries, depicted as the dotted lines in Figure 3.2. For example, the fully filled subcell Ω11 yields T11 = ∂ ΩW 11 ˆ + q · ndΓ = −q11 · ∂ ΩE11 ∂ ΩN 11 ˆ q · ndΓ ˆ + ndΓ (3.39) ∂ ΩS11 ˆ ndΓ where the superscript D ∈ {E, S,W, N} of ∂ ΩD i j specifies the subcell boundary. Note that in each subcell, the normals along the boundary segments are constant. Use Γ to represent the length of the boundary in the subcell (i, j). Equation 3.39 can lΩ i, j be simplified as T11 = −q11 · E l11 (3.40) S −l11 Similarly, for the partially filled cells, e.g. Ω22 , we get ˆ q · ndΓ T22 = Γc (3.41) = −q22 · ˆ + ndΓ Γa ˆ ndΓ , Γb Γ notation, Equation 3.41 where Γa , Γb , Γc are depicted in Figure 3.2. Use the lΩ i, j becomes T22 = −q22 · 19 Γa −l22 Γb l22 (3.42) Computing the other subcells in a similar manner and substituting v into the flux representation yields the linear constraint on the velocity: ∑ Ti j (v) = i, j 1 − Jn V. ∆tJn (3.43) Finally, each constraint of the form Equation 3.43 is assembled into a global system of constraints: GV v = c. (3.44) A full discretization of the momentum equation is given by M GVT GV v λV = pB c , (3.45) where M is a global mass matrix, pB are external impulses and λV are Lagrange multipliers. The above discretization is known to suffer from locking. This problem can be remedied via adding a stabilization term to Equation 3.45 to avoid spurious modes. We follow the approach by Misztal et al. [2012] which places εL in the (2, 2) block, where L is the Laplacian over the simulation mesh and ε is a small constant (see Section 5.3). 3.4 Contact and Collision Collision detection has been well studied especially in graphics. Note that our method does not require an explicit surface mesh for collision detection and contact. For example, Levin et al. [2011b] detects object interpenetration by subsampling. However, having a surface mesh explicitly enables the use of any existing collision detection engine, which could boost both performance and accuracy. As an example, we utilize a state-of-the-art collision detection library [Wang et al., 2012] in the implementation. As a bonus, the reconstructed surface mesh can directly be used for rendering purposes. 20 3.4.1 Surface Reconstruction We leverage our Eulerian mapping to perform surface reconstruction directly from the material coordinates X. This imparts a guarantee that the surface mesh will always return to its original shape when the mapping φ −1 becomes identical (something that cannot be said for advecting mesh vertices in an Eulerian flow). We store a surface mesh for our object in the material space. The Eulerian displacements stay small over the course of a time step. This allows us to perform a search for the position in the physical space x, for a given vertex of our mesh using the simple, iterative procedure given by Algorithm 2. Once this is done we can use the mesh for both collision detection and rendering. The availability of a fixed surface representation in X also allows us to perform texture mapping with zero drift. Algorithm 2 Search procedure for surface reconstruction. 1: 2: 3: 4: 5: 6: 7: 8: 9: for v ∈ Surface Mesh do p = vn e = v0 − X (p,tn+1 ) while ||e|| > ε do e = v0 − X (p,tn+1 ) p = p + Fe end while vn+1 = p end for In Algorithm 2 v0 is a vertex in our surface mesh in material space(usually, this is the initial configuration, vn and vn+1 are the vertices in the physical space at the time step tn and tn+1 , ε > 0 is a user defined tolerance and F is the Eulerian deformation gradient computed at p, which is the current guess of the solution. Note that at time tn+1 we do not know the physical space position of vn+1 so we guess that it has not moved since time tn . We compute a material space error using the function X(x,tn ) which is simply the Eulerian mapping and finally an physical space correction using F. Upon completion of the algorithm p holds the correct physical space position of v at tn+1 . In practice we execute Algorithm 2 in parallel for each mesh vertex. 21 3.4.2 Contact Simulation of multiple muscles in close contact requires formulating appropriate contact constraints, the details of which we now give. Consider a colliding pair of objects A and B with contact surface ΓAB . Non-interpenetration constraints are formulated on the velocity level so that the relative velocity in the normal direction is always nonnegative along the contact surface: vrel · nˆ ≥ 0. (3.46) In FVM we integrate Equation 3.46 in each grid cell Ωc where contact occurs to obtain ˆ vrel · ndΓ Γc ˆ ≥ 0, vA − vB · ndΓ = (3.47) Γc where Γc = Ωc ∩ ΓAB . By expressing velocities in terms of shape functions NA and NB Equation 3.47 becomes ∑ s = gTc v ˆ · vAs − ∑ NAs ndΓ s Γc ˆ · vBs NBs ndΓ (3.48) Γc ≥ 0. (3.49) Finally, the constraints are assembled into a global constraint Jacobian matrix GC = {gc }c=1,...,m (3.50) When muscles are in close contact, there exist connective tissues keeping the muscles from separating. Therefore we treat the inequalities in Equation 3.49 as equalities in the case of inter-muscular contact. This benefits the QP solver described in Section 5.3 as the number of inequality constraints is greatly reduced. As a bonus, this approach offers a simple and natural solution for simulating connective tissue without additional work using, e.g., the method by Teran et al. [2005]. 22 3.5 Plasticity One advantage of Eulerian simulations is that we can avoid material and physical space remeshing when simulating plastic deformations. Our implementation of elastoplasticity is based on that of Barteil et al. [2007] and Wicke et al. [2010]. The multiplicative plasticity model starts from the following decomposition of the total deformation gradient into elastic part and plastic part Ftotal = Fe F p , (3.51) where Fe is the elastic deformation gradient and F p is the plastic deformation gradient. The evolution of the plastic deformation gradient is given by tn+1 tn −1 (F−1 = (F−1 p ) p ) ∆F p , (3.52) −1 where ∆F−1 p is an update applied to F p . This update is computed using the singu- lar value decomposition (SVD) of the current elastic deformation gradient: Fe = UΣVT . (3.53) The plastic update is then computed as ∆F−1 p where γ = ν σ − σy σ Σ =V (detΣ)1/3 −γ VT , (3.54) , 0 ≤ γ ≤ 1. γ is a function of the flow rate ν and the plastic yield σy . Our approach is distinguished from previous methods [Bargteil et al., 2007, Wicke et al., 2010] by the manner in which the plastic deformation is stored. Though we allow Eulerian deformation, the availability of a mapping to X allows us to store the plastic deformation in the material space. When computing (Equation 3.51) in physical space, we perform a lookup for F p in X. When evolving F p , we construct a mapping from the material space to the physical space using the same method as surface reconstruction described in Section 3.4.1. The advantages of our approach are two-fold: a) it does not require remeshing and avoids artificial 23 plasticity; b) by storing F p in the material space, the plastic deformation gradient will not smear out during advection. Figure 3.3: Two plastic rods are colliding with a rigid rod. One example is shown in Figure 3.3. The two rods are moving towards the rigid rod at a high speed. After the collision, the ones under plastic deformation have their rest shape permanently changed and the contact area resembles the surface of the rigid rod. 24 Chapter 4 Eulerian-on-Lagrangian Simulation Eulerian simulations have proven to be indispensable tools in the computer graphics field. They are used to produce stunning animations of complicated fluids, and more recently, of solids and granular materials. Because these simulations rely on a fixed spatial discretization they avoid many of the perils encountered by Lagrangian simulation of highly deformable objects. Chapter 3 has given most of the fundamental details for simulating incompressible hyperelastic materials. However, these methods suffer from inherent drawbacks stemming from the fixed spatial discretization which is their defining characteristic. Chief amongst these is the necessary discretization of the entire simulation domain. For a musculoskeletal simulation, the muscles must be able to move together with the bones freely in space. Other difficulties include computing accurate and dissipation free advection, and potential time-step restrictions. Many algorithms attempt to avoid these issues by replacing parts, or all of the simulation with Lagrangian methods, but this is done at the peril of losing the advantages of the Eulerian approach. In this chapter, we present a clean, hybrid technique which allows us to leverage the benefits of Lagrangian and Eulerian simulators by treating rigid and low frequency modes of a deformable object as Lagrangian and the high frequency deformable modes as Eulerian. We will show that this approach alleviates or significantly reduces all of the problems listed above but still allows us to resolve large deformations conve25 niently in an Eulerian context. We also detail a contact-aware solver for these types of mixed Eulerian-Lagrangian simulations. For the reader’s convenience we will begin by outlining the notational conventions used in this chapter (Table 4.1). With a slight abuse of notation we denote both nodal variables for single elements as well as entire objects by the same bold face letters, and the meaning will be evident from the context. When necessary we denote the stacked vector of all configuration variables as q. Table 4.1: Notation used Notation f f • f f˙ [f] Definition scalar or 3-vector n × 1 vector resulting from stacking all f df dt df dt in spatial domain in intermediate domain cross product matrix: f× Material(z) Intermediate(y) Spatial(x ) Figure 4.1: The continuous domains used for the Eulerian-on-Lagrangian simulator as well as the discrete grid structure stores in the Eulerian domain. The mapping between the domains is described using a set of generalized configuration variables. Our exposition begins by describing three separate continuous domains: The material domain, the intermediate domain (also known as the reference domain) and the spatial domain (Figure 4.1). We define invertible mappings φi between these spaces and parameterize these mappings using generalized configuration vari26 ables qi . φ2 maps reference objects into their deformed state due to Eulerian motion and φ1 maps deformed objects into space via Lagrangian motion. The use of configuration variables is key to our approach since it allows us to build in any Lagrangian motion using a reduced set of coordinates. In practice we compute φ2−1 by advecting material coordinates on a uniform grid located in the intermediate domain. In order to elucidate the mechanics involved, we begin with the pedagogical case of a rigid Lagragangian motion. 4.1 Rigid Motion Eulerian-on-Lagrangian Simulation To begin we define mappings φ1 and φ2 as φ1 : x = R (y + p) (4.1) φ2 : y = z + u (z,t) , where R and p are a rotation and translation (in the intermediate space) that map y to x, and u is the displacement of a material point at z to its deformed position y. Thus φ1 represents the rigid motion (Lagrangian), while φ2 represents the deformation (Eulerian). By taking the time derivative of the composition of these functions we arrive at the spatial velocity of any point in the material domain • x = R [y]T ω + ν + v , (4.2) where [y] is the cross product matrix, ω is an angular velocity and ν is a linear velocity of the rigid motion and v = u˙ which results from Equation 4.1. In the intermediate space u is a function of y(t) and so u˙ corresponds to the material derivative, Du Dt . This gives rise to the Eulerian-Lagrangian velocity • x = R [y]T ω + ν + Du . Dt (4.3) We take the Lagrangian configuration variables, q1 , to be the axis-angle of rotation θ , i.e., R = exp[θ ], and position p (as seen in y) of the rigid frame in which the Eulerian simulation is embedded. The displacements are then discretized in the intermediate domain, y; this is 27 in contrast with the material domain in the Lagrangian approach and the spatial domain in the Eulerian approach. We use a hexahedral finite element grid with trilinear shape functions; the nodal displacements become our configuration variables q2 . For a single element this yields ω q˙ = ν . v θ q = p u (4.4) Again we note that the total time derivative of the Eulerian velocity implies the application of the material derivative. In our method this is computed using an independant advection step (as is usually done in fluid mechanics, see Bridson’s book [2008] for a good description). Other methods bake this material derivative into the equations of motion (see Sueda et al. [2011] for a depiction of such an approach). For the ith element we can define the Lagrangian function of our system as 1 L = q˙ T Mi q˙ −V (q) 2 where Mi = Ωi [y] [y]T ρ J¯ [y]T [y] I T NiT [y] NiT [y] Ni (4.5) Ni dΩi , NiT Ni (4.6) ρ is the density, J¯ is the deformation gradient determinant, Ni is the matrix of element shape functions and V is a potential energy accounting for all integrable forces. We evaluate these integrals numerically [Belytschko and Moran, New York, 2000]. It is convenient to think of this mass matrix in block form in which the diagonal blocks are the Lagrangian and Eulerian mass matrices while the off-diagonal blocks are coupling matrices (Figure 4.2). Note that, as in a reference coordinate formulation, the Eulerian displacements can be used to directly compute strain and avoid drift away from the object’s rest configuration. We pause here to offer a second derivation which demonstrates the general nature of the method. 28 4.2 Linear Modal Eulerian-on-Lagrangian Simulation In order to show how other Lagrangian discretizations can be utilized we will now repeat the process above for a more general case in which the Lagrangian motion is represented using linear modal models. We again begin by defining our mappings for an arbitrary mesh element, i: φ1 : x = N (y) Uq1 (4.7) φ2 : y = z + u (z,t) , where N is a matrix of shape functions, U is the matrix of linear modes for the element (comprised of the appropriate rows from the global linear deformation basis), and q1 are the degrees of freedom for the modal model. The spatial position of a given reference point is then x = N (z + u (z,t)) Uq1 (4.8) and the velocity can be computed as • x = F1 v + NUq˙1 (4.9) where F1 is the Lagrangian deformation gradient, given by ∇y qT1 UT NT and v is (as before) the intermediate space Eulerian velocity. In this case we choose our generalized coordinates to be q1 and the intermediate space, Eulerian displacements u. We form the Lagrangian of the system and write the per-element mass matrix, which in this case becomes Mi = 4.3 Ωi ρ J¯ UT NT NU UT NT F1 N NT FT1 NU NT FT1 FT1 N dΩi . (4.10) The Momentum Equation We note that both the rigid body and linear modal element mass matrices share a similar block structure. In fact this structure is common to all Eulerian-on- 29 Figure 4.2: The block structure of the Eulerian-on-Lagrangian mass matrix using rigid motion. The diagonal blocks are Lagrangian and Eulerian respectively while the off-diagonal blocks represent the coupling terms of the system. Lagrangian (EoL) derivations. The mass matrices appear as M1 M12 MT12 M2 (4.11) where the diagonal blocks are the mass matrices for our individual Lagrangian and Eulerian dynamical systems (respectively) and the off-diagonal blocks are the coupling matrices. Applying the Lagrange-d’Alembert principle [Lanczos, 1986] gives rise to the equations of motion for our element. After assembling the perelement mass matrices and force vectors we arrive at Mq¨ = fqvv + fe + fb + fc = f (4.12) where M is the global mass matrix, fqvv are the centrifugal and Coriolis forces acting on the Eulerian domain, fe are the forces of elasticity, fb are body forces, and fc are forces due to contact. The force fqvv results from the quadratic velocity terms in the Euler-Lagrange equations [Murray et al., 1994]. We observe that this global mass matrix retains the block structure of Equation 4.11 which is illustrated 30 in Figure 4.2. From here we can proceed in a completely general fashion, assuming only that the mass matrix has the described structure. ¨ which can be thought The crucial part of this system is the acceleration term q, of as the solution to the set of differential equations q¨1 ∂v ∂t + v · ∇v = M−1 f (4.13) where we have expanded the material derivative for the sake of clarity. One oftused method of solving such a differential equation is splitting, wherein one first solves q¨1 ∂v ∂t = M−1 f (4.14) and then uses the result, v∗ , to solve nodal equations of the form ∂ v∗ + v · ∇v∗ = 0 ∂t (4.15) which is now an advection and can be solved using any relevant algorithm (e.g., [Lentine et al., 2011]). Advection occurs on an Eulerian grid defined in the intermediate space. Algorithm 3 A high level overview of the Eulerian-on-Lagrangian solution procedure. Here we explicitly denote the partial time derivatives of Eulerian qualities for clarity. 1: Solve M q¨1 ∂v ∂t and v∗ = f for q˙ t+1 1 vt+1 = advect(v∗ ) q1 q˙1t+1 3: Solve ∂ u = for q1 and u∗ vt+1 ∂t 4: u = advect(u∗ ) 2: Thus our solution procedure is equivalent to a standard Lagrangian dynamics simulator with extra advection steps required at the velocity and position levels (Algorithm 3). Below we will describe our particular solver in greater detail. To solve step 1 of Algorithm 3, we discretize the acceleration q¨ and solve at 31 the velocity-impulse level, which gives def Mq˙ t+1 = ∆tf + Mq˙ t = pt+1 . (4.16) Using the block structure of M and suppressing the time step t + 1 for clarity, M1 q˙ 1 + M12 v∗ = p1 (4.17a) MT12 q˙ 1 + M2 v∗ = p2 , (4.17b) where M1 is the nonsingular Lagrangian mass matrix, M2 is the nonsingular Eulerian mass matrix, M12 is the coupling matrix between the two modes, p1 is the generalized Lagrangian momentum and p2 is the generalized Eulerian momentum(Figure 4.2). Because we lump the mass to the nodes, M2 is diagonal. 4.4 Exploiting Redundancies Solving Equation 4.16 is difficult because the mass matrix, M, contains a nontrivial nullspace. Consider the case of a simple translation. The same motion can be described by advection through the Eulerian grid or a Lagrangian translation of the grid itself (see Figure 4.3). The total velocity of the system is unique but its decomposition into Eulerian and Lagrangian parts is not. We can illustrate this more formally using the fundamental property of EoL: any momentum acting on the system can be solved purely on the Eulerian velocities v∗ . Thus, by setting q˙ 1 = 0 we get p1 = M12 v∗ p2 = M2 v∗ (4.18) ⇒ p1 = M12 M−1 2 p2 . Now suppose we induce a momentum via a Lagrangian velocity vL . Then p1 = M1 vL p2 = MT12 vL . 32 (4.19) Combining this with Equation 4.18 yields T M1 vL = M12 M−1 2 M12 vL (4.20) T ⇒ M1 = M12 M−1 2 M12 since vL was arbitrary. Finally, suppose we have a solution (q˙ 1 , v∗ ) to EquaT ˙ ). Evaluating the left hand tion 4.17b which we write as v∗ = M−1 1 2 (p2 − M12 q side of Equation 4.17a gives M1 q˙ 1 + M12 v∗ T ˙ 1 + M12 M−1 =(M1 − M12 M−1 2 M12 )q 2 p2 (4.21) =M12 M−1 2 p2 =p1 , which follows from Equation 4.18 and Equation 4.20. Thus a solution to Equation 4.17b always satisfies Equation 4.17a and it is therefore redundant. Examination of Equation 4.17b reveals an underdetermined system with a nullspace of dimension equal to the number of Lagrangian configuration variables. To see this, observe that any solution to Equation 4.17b may be written as q˙ 1 v∗ = I T −M−1 2 M12 q˙ 1 + Due to the full row rank of M12 , the matrix nullspace. In other words, we are free to choose 0 M−1 2 p2 . I T −M−1 2 M12 any q˙ 1 and the (4.22) is a basis for the Eulerian velocities “compensate” via Equation 4.22. We specify a unique solution by maximizing the size of the time step taken which is equivalent to minimizing ||v∗ ||. Since the time step is specified by the maximum grid velocity, the ∞-norm would be ideal; however, this requires the solution of a linear program on the order of the number of Eulerian variables. Instead 33 Figure 4.3: Due to redundancy of the Lagrangian modes, the motion of an object can be specified by advection through the fixed Eulerian grid (top right), or fixing the object in the grid and moving the grid itself via the Lagrangian degrees of freedom (bottom right). we approximate with the 2-norm, giving −1 −1 T ˙ 1 ||2 , q˙ t+1 1 = argmin||M2 p2 − M2 M12 q (4.23) −1 T t+1 ˙1 . v∗ = M−1 2 p2 − M2 M12 q (4.24) q˙ 1 This is a least squares problem (LS) in the small number of Lagrangian variables, which can be computed efficiently. As M2 is diagonal its inverse is computed explicitly. In this discussion we have assumed an Eulerian momentum p2 is given. In Section 4.5 we detail a method for computing the global Eulerian velocity v in the spatial domain. We then simply compute p2 = M2 Ov, where O selects the components of v for the given object and transforms its spatial velocity to its intermediate velocity, since this is the space in which advection occurs. Finally, it will be useful for us to rewrite Equation 4.24 as v∗ = projMC Ov, T where MC = M−1 2 M12 and proj is the projection operator. 34 (4.25) 4.5 Solving the Dynamics in the Presence of Contact We first solve the dynamics exclusively on the Eulerian grid (i.e. only using the Eulerian degrees of freedom), which can be done using any Eulerian simulation method suitable for resolving contact. In our case, we invoke Gauss’ principle of least constraint [Lanczos, 1986], which gives the constrained optimization problem 1 T ˜ a Ma − ˜fT a 2 subject to a ∈ A , minimize (4.26) ˜ is a global mass matrix, a and ˜f are the accelerations and forces defined where M in the spatial domain, and A is a constraint manifold imposed on the acceleration; these constraints ensure interpenetration between objects does not occur. Notice ˜ is assembled from the Eulerian mass matrices M2 in Equation 4.17. First that M order discretization of the acceleration gives rise to a quadratic program (QP) on the spatial velocities with constraints we now derive. Consider a colliding pair of objects A and B with surfaces ΓA and ΓB , respectively. Due to the discrete time integration there is some volume of intersection Ω. We impose a velocity level constraint specifying that the volume of intersection at the next time step must be smaller than Ω, giving the surface integration constraint vA · ndΓ vB · ndΓ ≤ 0, + ΓA ∩Ω (4.27) ΓB ∩Ω where n are surface normals. Spatial velocities are represented on the Eulerian grid by trilinear shape functions. Thus expanding the velocity terms in Equation 4.27 yields the constraint n ∑ vAs s=1 NAs · ndΓ ΓA ∩Ω + vBs NBs · ndΓ ΓB ∩Ω T =j v≤0 (4.28) where s indexes the nodal velocities and their associated interpolating functions given by NAs and NBs . 35 The constraints are assembled into a global constraint matrix GC and the QP solved at each time step is given by 1 T ˜ v Mv − p˜ T v 2 subject to GC v ≤ 0, minimize (4.29) where p˜ is a globally assembled momentum vector. In the case of rigid bodies with externally prescribed motion the formulation changes slightly as the right hand side is no longer 0. Collision detection and the computation of the surface integrals in Equation 4.28 require a representation of each object’s surface in the spatial domain. Any suitable method could be used, such as the volumetric method of Levin et al. [2011b] or recent methods based on ray tracing [Wang et al., 2012]. In our implementation we use the reconstructed surface (see Section 3.4.1), with a ray tracing collision detector to produce a collection of intersection rays, barycentric coordinates w.r.t. the surface mesh, and normals for every pair of colliding objects. Since spatial velocities are stored in the intermediate domain (on the Eulerian grid), we use the barycentric coordinates of the intersected rays to find the corresponding positions y and their associated degrees of freedom, from which the shape function values can be computed as well as the nonzero pattern of J. Finally, we note that constraints between objects are decomposed via a collision grid defined on the spatial domain to enhance the resolution of constraints along deforming surfaces. Equation 4.29 is solved via an efficient primal-dual active set method [Ito and Kunisch, 2008]. Due to the diagonal structure of the mass matrix, the complexity of the QP solver is O(m3 ), where m is the number of constraints. Once the spatial velocity v is known, we can compute the Eulerian momentum p2 of each object. Thus step 1 of Algorithm 3 amounts to a QP solve to determine the global momentum followed by a LS solve to optimally compute Eulerian and Lagrangian velocities. See Table 6.1 for computational runtimes. Although the method presented here handles most of the contact scenarios well, there are certain limitations. Due to volumetric representation of objects, handling the contact between solids and thin shells is non-trivial in the current simulation framework. See [Faure et al., 2007] for a generic and efficient method based on 36 an Eulerian contact model and [Robinson-Mosher et al., 2008] for the coupling between Eulerian fluid and shells. The other limitation is the difficulty of resolving contact at sub-cell level. For example, when a needle is in contact with a soft object, our method is not able to capture all the details around the needle tip. 37 Chapter 5 Musculoskeletal Simulation and Activation Figure 5.1: Muscle, tendon and bone structure The musculoskeletal system is an extremely complex system comprised of muscles, bones and tendons. In our simulation, we treat bones as rigid bodies, muscles as actuators attached to bones, and the tendons as inextensible material, kinematically prescribed by the bones. To elucidate our approach further, consider the human arm example as depicted in Figure 5.1, consisting of the biceps attached to the Humerus. Arm kinematics are characterized in biomechanics using the concept of a “moment arm;” conceptually this is equivalent to the radius of a pulley approximating the joint structures on which the muscle’s tendon moves (in this case, the elbow). The biceps and forearm muscles are connected via tendons. Tendons are much stiffer than muscles and can be well approximated as inextensible; 38 this approximation is implicit in biomechanical simulators such as OpenSim [Delp et al., 2007]. Hence their kinematics are determined entirely by bone configuration. Since the elbow is treated as a revolute joint, tendon excursion is simply d = θ r, where θ is the rotational angle and r is the moment arm. Finally, the biceps muscle serves as the actuator for the lower arm, which we model as two separate muscles since they have different origins and hence different mechanical actions. Even though there is some lateral force transmission between adjacent muscles, this appears to be small in the physiological range [Maas and Sandercock, 2010]. Data such as moment arms of each tendon are taken from [Holzbaur et al., 2005]. We now describe the approach that we follow for simulating bones, tendons and attachment constraints. Finally we explain how to solve the global optimization problem at each time step. 5.1 Bones and Tendons Bones are typically modeled as rigid bodies. Joints connecting those bones ensure smooth relative motion and are critical for the musculoskeletal simulation. Methods modeling the joints in an anatomically correct way widely exist in biomechanics. In graphics, Lee et al. [2008] has developed the spline joint algorithm for exploring the underling configuration space for one DOF joint such as the knee. Although it presents an interesting mathematical model, the method needs further biomechanical validation. Its extension to multi-DOFs joints is not straightforward. Alternatively in practice, combinations of basic joint types are heavily used in robotics and graphics. In this thesis, we build the skeletal model based on basic joint types and reduce the bone dynamics to a problem of solving a multi-body system. We use maximum coordinates for the bone dynamics [Cline, 2002]. Denote q˙ = (x, ˙ y, ˙ z˙, α, β , γ) ∈ R6 as the rigid body velocities. The first three components (x, ˙ y, ˙ z˙) represent linear translational velocities and the other three (α, β , γ) represent angular velocities. Those velocities are solved by invoking Gauss’s principle of least constraints as described in Section 3.2. After discretizing the joint accelerations using forward 39 Euler, we arrive at the following QP: 1 T q˙ Mq˙ − q˙ T h 2 q˙ n+1 =argmin q˙ (5.1) subject to GB q˙ = b where h is the predictive momentum of the rigid bodies. Joints are modeled by adding equality constraints GB to Equation 5.1. Details of forming the constraints for revolute joints, ball and socket joints and others can be found in Cline’s thesis [2002]. Denote ln+1 = (xn+1 , yn+1 , zn+1 )T as the linear translation at the time step tn+1 . Once the velocities are solved from Equation 5.1, the DOFs are updated using forward Euler ln+1 = ln + dt ˙l Rn+1 = Rn ∆R (5.2) where ∆R is a small rotational matrix caused by the angular velocities (α, β , γ) over the time step ∆t. The rigid body transformation matrix Tn is assembled as Tn+1 = Rn+1 ln+1 0 1 (5.3) Our tendon model is based on the assumption that tendons are inextensible since they are very stiff comparing to the muscles. This assumption bears the merit that it is both biomechanically sound and computationally efficient. The dynamics of the tendons can be simulated by one dimensional models such as strand [Sueda et al., 2008]. It is also common to model the tendons as quasi-static, so that the full motion can simply be described in terms of moment arms. We illustrate this idea by showing an example of a rigid body system with two revolute joints. In Figure 5.2, the two links are initially placed horizontally (the bottom row). Later they are rotated by θ1 and θ2 . The tendon is denoted by the thick and orange line. Its one end is inserted into the second link and the whole tendon lies on the surface of the rigid links. Denote r1 and r2 as the radii of the two joints (also the 40 Figure 5.2: A 2D rigid body system with an inextensible tendon (orange line) moment arms). The excursion e at the far end of the tendon is e = θ1 r1 + θ2 r2 = r1 r2 θ1 (5.4) θ2 In short, we write Equation 5.4 as e = Lθ , (5.5) where θ is the stacked variable for the rotation angles (θ1 , θ2 )T . The linear velocity at the far end of the tendon is e˙ = Lθ˙ , (5.6) If an external force f is applied to the tendon, the corresponding torques ap- 41 plied to the joints are τ = LT f (5.7) In Section 5.2, we use the results from Equation 5.5 and Equation 5.6 to glue the muscles to the tendons. 5.2 Coupling of Bones, Tendons and Muscles First, we couple muscles with the bones and tendons via attachment constraints. Typically, the proximal (origin) end of a muscle attaches over a broad area of a bone, and the distal (insertion) end has a more prominent tendon. The Brachialis, for example, has a large area of origin on the Humerus, and the distal end connects to the forearm via tendons. For a muscle with attachment boundary Γ either at the origin or at the insertion, we enforce the Eulerian velocity on Γ to be the same as the tendon velocity e˙ , which is the relative velocity to the EoL domain used to describe the muscle. vn+1 (x) = e˙ x ∈ S, (5.8) where S ⊆ Γ is the attachment area. Its position in physical space is determined by excursion e and rigid body transformation. Writing v in terms of shape functions Ni and integrating Equation 5.8 over a grid cell Ωc yields ∂ Γ = e˙ ∑ Ni vn+1 i Γ∩Ωc i ∂ Γ, (5.9) Γ∩Ωc Note that in the case of muscle-bone attachments, e˙ on the attachment area is zero throughout the simulation as the bones are fixed with respect to the simulation grid. A muscle normally has two attachments, one for origin and the other for insertion. We globally assemble these constraints, giving GA v = a, (5.10) where GA is the attachment constraint matrix and a represents the right hand side. We modularize the musculoskeletal system to separate the tendon-skeleton subsystem from the muscle subsystem. This provides freedom to “plug-in” submodule solvers seamlessly, e.g., more sophisticated rigid body dynamics engines 42 than the one described here. We enforce weak coupling between the subsystems to achieve this modularization. We detail the systems below. Equation 5.1 has given the QP for the bone dynamics. Here with tendon forces coupled, the new QP for skeleton-tendon is. q˙ n+1 =argmin q˙ 1 T q˙ MS q˙ − q˙ T (pS + pCn ) 2 (5.11) subject to GS q˙ ≤ d, where MS is the inertia tensor, pS is the external impulse and pC is the coupling impulse resulting from attachment constraints. Note that all quantities defined are globally assembled. The constraints of Equation 5.11 are general. The inequality constraints can enforce, for example, joint angle bounds. Similarly, the muscle QP is vn+1 =argmin v 1 T v Mv − vT pB 2 (5.12) subject to GV v = c, GC v ≥ 0, GA v = a, where pB is the impulse due to advection and all the forces. M is the global mass matrix. We simplify the numerics by lumping mass to the nodes so that M is diagonal. After Equation 5.11 is solved, we get q˙ n+1 and subsequently formulate the attachment constraints using the newest rigid body DOFs and Equation 5.12 can be solved. Note that the coupling impulses λAn+1 have been implicitly computed as the Lagrange multipliers associated with the attachment constraints: pCn+1 = GTA λAn+1 . (5.13) Once we know pCn+1 , we can advance to tn+2 by solving Equation 5.11 using pCn+1 . A high level overview of our solution procedure for a single time step is given in Algorithm 4. 43 Algorithm 4 Solution Procedure 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 5.3 Initialize λA0 = 0, q0 = 0, v0 = 0, t0 = 0 and n = 0 while tn < ttotal do Compute pCn = GTA λAn Solve the skeletal QP for q˙ n+1 Update qn+1 = qn + ∆t q˙ n+1 Advect vn to obtain v∗ Formulate the attachment constraints in Equation 5.12 using qn+1 Solve the muscle QP for vn+1 and λAn+1 Xn+1 = advect(Xn ) tn+1 = tn + ∆t n = n+1 end while Solver Each time step requires the solution of the skeletal and muscle QPs, given by Equation 5.11 and Equation 5.12, respectively. Both are solved via a primal-dual activeset method [Ito and Kunisch, 2008]. Each iteration of the active-set method requires the solution of the following symmetric, indefinite linear system or KKT system. For the skeletal QP, the size is small and the solution is computed via a direct solver. The KKT system of the muscle QP is M GVT ˜T G C GTA GV −εL G ˜C GA pB λV c = λ 0 C λA a v∗ (5.14) or Ku = b, ˜ C consists of rows from the contact Jacobian GC corresponding to active where G constraints. Due to the size of Equation 5.14, we employ MINRES, a Krylov subspace iterative method for solving large sparse symmetric linear systems. Iterative methods only require matrix vector products of the form Kx, avoiding the need to compute the matrices in Equation 5.14 explicitly. Crucial to the efficiency of Krylov subspace methods is the choice of a preconditioner matrix P, which solves 44 the equivalent system P−1 Ku = P−1 b. (5.15) MINRES requires a positive definite preconditioner, and we use the SPD matrix M P= L , S (5.16) where S= ˜ T M−1 G ˜C G C GTA M−1 GA . (5.17) From (5.15), the matrix vector products required in MINRES requires a solve with P, which amounts to solving with the matrices M, L and S. Solving Mx = b with the diagonal mass matrix M is trivially computed in place on the grid of each object. As the contact and attachment constraints occur on the surfaces of objects, ˜ C and GA explicitly, S is expected to remain relatively small. Therefore we form G from which S can be efficiently computed. It can be shown that S exhibits a banded structure, and we solve systems of the form Sx = b via a banded Cholesky factorization S = RRT . Finally, we never form the matrices GV and L explicitly; rather, we solve systems of the form Lx = a (5.18) via an inner iterative solver: the conjugate gradient (CG) algorithm. CG, another Krylov subspace iterative method, typically works well on Laplace’s equation, but it is a nonstationary solver. Using CG for parts of the preconditioner violates the convergence theory for MINRES. In practice, the presented inner/outer iteration scheme converges to the solution, but the number of inner CG iterations barely affects the outer MINRES iterations. We found that a simpler preconditioner M P= I , S 45 (5.19) exhibited better performance. For future work, we are going to use flexible Krylov methods such as FGMRES to iteratively solve Equation 5.14 with CG as the preconditioner for the L block. Table 6.2 shows the performance of our solver for various problem data encountered in a typical simulation. All factorizations and solves are performed using the GPU package CULA [Humphrey et al., 2010]. 5.4 Data-driven Muscle Activation The essential idea of the muscle activation model is the change of reference configuration. We follow the same approach as simulating elasto-plastic materials (Section 3.5) by decomposing the total deformation gradient into elastic part and plastic part. Ftotal = Fe F p , (5.20) ˜ in material space. Using For muscle activatoin, we store the active configuration X Fa for the active deformation gradient, we treat it like F p as in plasticity simulation and write Fe as Fe = Ftotal F−1 a (5.21) Figure 5.3: The mapping between active space, material space and physical space ˜ X and x through the Figure 5.3 shows the relationship between the spaces X, −1 ˜ = ψ(X) = ψ(φ (x)). mapping X In our method, it’s unnecessary to compute Fa explicitly. Instead, we directly 46 ˜ by the mappings ψ and φ −1 and then compute Fe as interpolate X Fe = ∂x = ˜ ∂X ˜ ∂X ∂x −1 (5.22) For Eulerian methods, this approach avoids tensor field interpolation or advection on the grid, thus results in a simpler algorithm and better results. We didn’t observe active deformations being smeared out or any artificial deformation due to activation. These benefits are brought about by the fact that, a priori, we know ˜ exists and is provided as input. Meanwhile we assume no residual stress in the X active configuration. The method shares some common elements with the work by Wicke et al. [2010], where they try to store as much plasticity as possible in material space by remeshing so as to minimize artificial plasticity. Since they are simulating elasto-plastic objects, they have to deal with residual stresses. One feature of this method is that we give modelers the option of choosing desired active shape freely. They could either manually deform it, as in Figure 5.4, or they could leverage any existing biomedical data. For example, the human ocular muscle activation simulated by our methods has successfully used MRI reconstructed data from the work by Wei et al. [2008]. Here, we do not impose any requirements on the volume of the activated configuration. As a side effect, the method could simulate how muscle remodels itself such that chronic muscle deformations do not generate stress [Bottinelli and Reggiani, 2006]. As in graphics, this provides a handy tool for deforming muscles exaggeratedly within a small amount of time. Figure 5.4: Muscle deformation by manually specified activation shape 47 During the activation process, due to a sudden change of the reference shape, the simulation requires tiny time steps. In muscle mechanics, such sudden activation implies adopting a step function in time, which is unrealistic. To enable smooth transition from relaxation to activation and to allow larger time steps, our method could take a one-dimensional function G and approximate intermediate ˜ − X). activation configuration as X + G (t)(X 48 Chapter 6 Results and Examples All the examples were simulated on an Intel i7 2.8 GHz CPU with an Nvidia GeForce GTX 580 graphics card. 6.1 Volume-Preservation To validate the effectiveness of the volume-preserving constraint with the presented stabilization method. We simulated crushing an extremely soft cube in 3D space Figure 6.1. The cube is placed in a huge horizontal force field. The force magnitude is uniform over the whole space, but the direction of the force depends on the x coordinate. More specificly, the force field is Ce , x1 < 0 1 f= ,C > 0 −Ce , x >= 0 1 1 (6.1) where C is a positive scaler and e1 = (1, 0, 0)T is a unit vector. The volume change of the objects is depicted in Figure 6.1. The three lines correspond to three simulations under investigation. Obviously the non-volume-preserving object (the red line) loses volume from the very beginning and is quickly smashed to a sheet. In contrast, volume-preserving solids (the green and blue lines) are more resistant against volume change under the external force field. The ∇·v = 0 formulation (the green line) loses volume over time, which is obvious in Figure 6.1. Our method with the proposed positional stabilization (the blue line) exhibits nice behavior. Its 49 Figure 6.1: The figure shows the volume change for three simulations in comparison. y-axis: object volume in percentage; x-axis: total simulation time in seconds. The red line is for the non-volume-preserving solid while the green and blue lines are for volume preserving solids. The green one is based on ∇ · v = 0 condition and the blue one is based on our method volume stays almost constant and the maximum volume change does not exceed 0.2% of the original volume. This result shows that our method does not accumulate numerical drift over time. Figure 6.2 shows that, without positional stabilization, once the force field is released, the volume-preserving solids do not necessarily return to their original shapes. As a reference, the blue line represents our method. The other three lines represent objects with different stiffness. The force field is applied from the beginning. After 1.2 seconds, it’s released. The objects start popping back to their original shapes, but the ones without positional stabilization (green, red and pur50 Figure 6.2: The figure shows the volume change for objects with different stiffness. The blue line representing our method is for reference. The object’s stifness is 1000kP. The other three lines represent objects with Young’s Modulus 1000kP, 1250kP and 1500kP respectively, simulated without stabilization ple) are experiencing difficulty during the recovery. In the examples, we have two important observations. First, at the end of the simulation, all three objects fail to return to the original volumes. Second, their volume losses do not correlate with the stiffness. For example, the stiffest object (the red line) eventually gains more volume than the softest object but less volume than the medium one. Based on these observations, we conclude that the presented positional stabilization is an effective method for preserving the actual volume. Without it, velocity-based volume-preserving methods tend to produce unexpected results and numerical drift easily plagues the simulations. Figure 6.3 demonstrate that contact and volume-preservation constraints are simultaneously resolved. One volume-preserving bunny (left) and a purely elastic bunny (right) are dropped to the ground at the same speed. Figure 6.3 captures 51 Figure 6.3: Comparison of volume-preserving and non-volume-preserving bunnies dropping to the ground. Left: volume-preserving bunny. Right: purely elastic bunny the moment when the two bunnies touch the ground and undergo large deformation. It’s visually clear that the non-volume preserving bunny has experienced a significant volume loss comparing to the volume preserving one. 6.2 Eulerian-on-Lagrangian Simulation Figure 6.5 shows the benefits of Eulerian-on-Lagrangian simulation for simple translation. Note that the Eulerian case (top row) requires discretization of all of space. This has two obvious limitations: a drastic increase in memory use and the restriction of the object’s motion to the domain. Note that the EoL simulation is free to roam wherever it pleases while wasting far less memory on extraneous grid points. The numerical advantages are made clear during rotation. When a constant angular velocity is applied to a purely Eulerian object it can be damped out by elastic forces on the grid (Batty et al. [Batty and Bridson, 2008] shows this for viscous fluids). When the rotational DOF is Lagrangian, damping is greatly reduced. Given an initial angular velocity the purely Eulerian case comes to rest while the EoL method continues turning almost indefinitely (Figure 6.6). Our solver automatically extracts these Lagrangian modes from arbitrary motion. In some sense the EoL method serves as a generalization of floating frame-based methods for purely Lagrangian simulation [Barbiˇc and Zhao, 2011, Galoppo et al., 2007] wherein a deformable Lagrangian simulation is coupled with a rigid frame. Explicitly exploiting the singularities in the EoL mass matrix allows us to utilize 52 any kinematic representation for the Lagrangian component such as linear modal models (Figure 6.4) Figure 6.7 shows our dynamic time stepping scheme in action, which ensures stability for the latter advection scheme. Once objects exit the collision they gradually become Lagrangian. At this point the time step returns to its maximum value which is determined by the (typically less restrictive) stability conditions of the Lagrangian time integrator. The dynamic time stepping and the continuous transition between Eulerian and Lagrangian states are handled automatically by the EoL solver. Figure 6.8 shows that the Eulerian-on-Lagrangian simulator enjoys a large advantage in terms of time step size over its purely Eulerian counterpart. This manifests itself as an order-of-magnitude increase for most of the simulation. Note that the dissipative effect of Eulerian advection, in the purely Eulerian case, can also be seen; the collision begins later in time, indicating that the object has lost some velocity. The method is well suited for simulating highly deformable objects in an unbounded domain. Furthermore, like its Eulerian relatives it does not require an object specific discretization thus avoiding the complexities of mesh generation. Timings are reported in Table 6.1. In our current implementation each object is processed serially, so the runtime for the least squares solver is proportional to the number of objects. For complex scenes, more collisions yield greater number of constraints and the QP solver performs accordingly. Memory requirements are largely due to the use of the ray-tracing based collision detection engine [Wang et al., 2012]. 53 Example 2 Spheres 2 Penguins Spinning Bar Plastic Cylinders Poking 16 Penguins Sphere Squeezer Toys Squeezer Dim 2 × 32 × 32 × 32 2 × 32 × 32 × 32 64 × 32 × 32 2 × 64 × 32 × 32 128 × 64 × 64 16 × 32 × 32 × 32 8 × 48 × 48 × 48 8 × 48 × 48 × 48 Lag DOF 6 6 16 12 12 6 9 9 Mem(MB) 799 777 737 833 931 1070 1018 1112 Steps 2773 1531 2087 1636 2284 4382 1918 2881 Forces 2.0 1.7 0.7 2.7 3.3 25.4 14.9 13.8 Collision 15.7 169 109.0 102.4 103.4 894.9 163.3 424.2 Assembly 2.9 1.9 0.8 2.4 64.1 168.5 166.6 207.9 LS 21.6 15.0 6.1 18.5 45.0 150.1 89.7 55.4 QP 54.4 22.3 12.0 28.4 27.4 93.4 538.6 159.4 Table 6.1: For each example, the size of the spatial grid (Dim), maximum GPU memory used, total number of time steps taken, and average runtimes in milliseconds for each stage of the algorithm is shown. 54 Dimensions 6 × 32 × 32 × 32 6 × 64 × 64 × 64 6 × 128 × 128 × 128 DOFs 3546 13851 68529 Volume 533 2750 16624 Attachment 147 279 627 Contact 38 111 181 MINRES Iters 68 62 56 LO (ms) 240 347 629 PC (ms) 744 1501 3969 total (ms) 279 484 1589 Table 6.2: Statistics for a typical time step of the iterative QP solver for varying dimensions of the human arm example. For each resolution, The number of DOFs and constraints of each type is shown, as well as runtimes of applying the KKT linear operator and preconditioner, in milliseconds. Figure 6.4: Left: Purely Eulerian motion. Middle: Purely Lagrangian motion. Right: Overall motion Figure 6.5: Top Row: Purely Eulerian translation. Bottom Row: Eulerianon-Lagrangian Translation 6.3 6.3.1 Musculoskeletons Eye muscles As shown in Figure 6.9, we simulate the extraocular muscles and the eyeball using the musculoskeletal simulator. The left pictures in Figure 6.9 are the raw MRI ˜ Wei et al. [2008] has images. They also specify the active configurations X. developed algorithms to extract muscle shapes from these images. The middle picture in Figure 6.9 is the eye that we model. The muscles are in their reference configurations in straight-ahead gaze, which we get from [Wei and Pai, 2008]. 55 E−on−L Eulerian Angular Velocity (m/s) 3 2 1 0 0 2 4 6 8 Time (s) 10 12 Figure 6.6: The angular velocity of an Eulerian-on-Lagrangian object (E-onL) and a purely Eulerian object undergoing rotation. Both were given an initial angular velocity of π radians/s. Step size 0.0167 0 0 50 100 150 Iteration 200 250 300 350 Figure 6.7: The time step, in seconds, of an Eulerian-on-Lagrangian simulation. Note that during contact the time step is automatically reduced. After contact the time step increases until it reaches the rendering frame rate. Among the six extraocular muscles, we only simulate Lateral Rectus (LR) and Medial Rectus (MR) and model the eyeball as a rigid body with one rotational DOF. Then we activate LR and MR one at a time. The right side of the figure shows the results. When the simulation reaches equilibrium, the activated muscle shapes on the right resemble the real data on the left except that the rotational angle is slightly smaller. One possible reason is that we do not model all the tissues for 56 E−on−L Eulerian Time Step (s) 0.0167 0.0016 0 0 1 2 3 4 Time (s) 5 6 7 8 9 Figure 6.8: The timestep size for both the purely Eulerian and Eulerian-onLagrangian simulations for the two ball collision. Note that in most cases the EoL time step is an order of magnitude greater than that of the purely Eulerian simulation. the eye. The other possibility is that the active shape from MRI is under external constraints, so that it does not truly reflect the active configuration and results in smaller forces. One quick remedy for this is extrapolating the active shape linearly. More biomechanically grounded methods are worth future investigation. 6.3.2 Arm Figure 6.10 shows the simulation of arm flexing. In essence, the arm that we model is a three DOF rigid body system actuated by six individual muscles. We use six different soft objects to simulate the muscles on the Humerus, two objects for the Biceps, three for the Triceps and one for the Brachialis. Each object has its own origin and insertion, so that they behave like separate muscles. As for the bones, we use four rigid bodies to model the arm structure. The Ulna and the Radius for the forearm are treated as one rigid body, so are the small bones in the hand. This simple arm structure is composed of three revolute joints representing the shoulder, the elbow and the wrist. The active shapes of the muscles in this example are manually drawn to accommodate artistic needs. Once the Biceps and the Brachialis fully activate, the arm flexes. We can clearly see from Figure 6.10 that the acti- 57 ˜ middle:the refFigure 6.9: Left: the MRI data of the active configuration X; erence configuration X; right: simulations of the activated muscles in physical space x vated muscles are deforming towards their active shapes. The Brachialis shortens and the Biceps bulges due to volume-preservation. As for the Triceps, they are sliding on the Humerus because the bone configuration is changing but their overall shapes remain the same since they are not being activated. Noting that significant improvements and extensions may be added, the current simulator has fulfilled the most basic functionality of a musculoskeletal system. 58 Figure 6.10: Up: the upper arm in the rest shape; Down: the Biceps and the Brachialis are activated 59 Chapter 7 Conclusions and Future Work 7.1 Limitations and Future Work Our Eulerian integrator treats elastic forces explicitly. Though dynamic time stepping provides stability, performance can suffer when simulating very stiff materials. Implicit integrators for these types of dynamic algorithms would be a useful avenue of work. Other fascinating areas to explore include simulating fluid-solid interaction in a unified Eulerian way, or simulating large fluid flows using many simultaneous Eulerian-on-Lagrangian domains. We also found it difficult to scale the musculoskeletal simulation. First, respecting the inherent complexity of musculoskeletons, the resulting system to be solved is highly constrained. Even the tailored QP solver (Section 5.3) consumes a considerable portion of the total running time. Taking the cost of collision detection and sub-sampling into account, we need to improve the performance for larger scale simulations. Second, modeling musculoskeletons requires more expert knowledge of anatomy, which is a hurdle for computer scientists. Fortunately, this can easily be overcome by collaborative efforts. Despite the fact that the simulator does not easily scale, it still opens up fruitful future directions. A direct extension is to simulate the full arm, including the muscles on the shoulder and on the forearm. The method can also be applied to other parts of the human body such as the leg, the torso and even the full body. It’s also interesting to couple the current system with a biomechanically based tendon 60 network simulator and to animate realistic hand motion by muscle activation. Incorporating frictional contact into the Eulerian solids method is another exciting future direction. Once equipped with friction, our simulator can be used to simulate the fingertip, which enables further constructive research in touch sensation and neural control. 7.2 Conclusion This thesis presents solid simulation techniques and its applications in musculoskeletal modeling and simulation. In terms of numerical methods, the major contributions are twofold. First we give a general formulation for Eulerian-onLagrangian simulation which combines the advantages of both Eulerian and Lagrangian discretizations. The introduced Lagrangian DOFs make Eulerian solids simulation more practical because the Lagrangian modes capture the global motion of the object and make the remaining Eulerian deformation small. This method brings the advantages of large time steps, small dissipation, and unbounded simulation domains to existing Eulerian methods. The other contribution is extending the Eulerian solids method to elasto-plastic material and volume-preserving material. Since the Eulerian solids method does not require remeshing in physical space, it avoids artificial plasticity. Moreover, by storing the plastic deformation gradient in material space, plastic deformation won’t smear out during long simulations. For volume-preserving solids, we take a fluid mechanics approach by enforcing the incompressibility constraint on the velocity-level. The presented stabilization method for tracking numerical drift is simple and effective. The examples show that object volume is well preserved both locally and globally over time. The second part of the thesis is devoted to the application of the developed methods to musculoskeletal simulations. Muscle simulation is based on the Eulerianon-Lagrangian method. Its rigid motion in space is specified by the bone dynamics while large deformation and close contact are robustly handled by the Eulerian method. By a few reasonable simplifications, the musculoskeletal simulator fully couples the bones, tendons and muscles. A novel muscle activation model that features easy implementation is presented. It follows the multiplicative decomposition law which was originally used for elasto-plastic simulation. We change 61 muscles’ material configuration during activation and this has the advantage that existing biomedical data such as ultra-sound and MRI can be used for the active configuration. Lastly, we have demonstrated an example based on real MRI data for the eye muscle simulation. Overall, the thesis has made several useful contributions to solids simulation with applications to musculoskeletal simulation and muscle activation. 62 Bibliography J. W. Banks, D. W. Schwendeman, A. K. Kapila, and W. D. Henshaw. A high-resolution Godunov method for compressible multi-material flow on overlapping grids. J. Comput. Phys., 223(1):262–297, 2007. ISSN 0021-9991. → pages 4 J. Barbiˇc and Y. Zhao. Real-time large-deformation substructuring. In ACM Transactions on Graphics (TOG), volume 30, page 91. ACM, 2011. → pages 52 Adam W. Bargteil, Chris Wojtan, Jessica K. Hodgins, and Greg Turk. A Finite Element Method for Animating Large Viscoplastic Flow. ACM Trans. Graph., 26(3):16:1–16:8, July 2007. → pages 6, 7, 23 Philip T. Barton and Dimitris Drikakis. An Eulerian method for multi-component problems in non-linear elasticity with sliding interfaces. J. Comput. Phys., 229 (15):5518–5540, 2010. ISSN 0021-9991. → pages 4 Christopher Batty and Robert Bridson. Accurate viscous free surfaces for buckling, coiling, and rotating liquids. In Proceedings of the 2008 ACM/Eurographics Symposium on Computer Animation, pages 219–228, July 2008. → pages 52 F Ben Belgacem, P Hild, and P Laborde. The mortar finite element method for contact problems. Mathematical and computer modelling, 28(4):263–271, 1998. → pages 8 T.B. Belytschko and J.M. Kennedy. Computer models for subassembly 63 simulation. Nuclear Engineering and Design, 49(1-2):17 – 38, 1978. ISSN 0029-5493. doi:10.1016/0029-5493(78)90049-3. URL http://www.sciencedirect.com/science/article/pii/0029549378900493. → pages 6 Ted Belytschko, Liu Wing Kam, and Brian Moran. Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, 2000. ISBN 471-98773-5. → pages 6 W. K. Belytschko and B. Moran. Nonlinear Finite Elements for Continua and Structures. Wiley, New York, 2000. → pages 28 Silvia S. Blemker and Scott L. Delp. Three-dimensional representation of complex muscle architectures and geometries. Annals of Biomedical Engineering, 33(5):661–673, May 2005. → pages 6 Javier Bonet and Richard Dixon Wood. Nonlinear continuum mechanics for finite element analysis. Cambridge university press, 1997. → pages 15, 16 R. Bottinelli and C. Reggiani. Skeletal muscle plasticity in health and disease: from genes to whole muscle. Springer, 2006. → pages 47 Robert Bridson. Fluid Simulation. A. K. Peters, Ltd., Natick, MA, USA, 2008. ISBN 1568813260, 9781568813264. → pages 11, 28 Mark Carlson, Peter J. Mucha, R. Brooks Van Horn, III, and Greg Turk. Melting and flowing. In Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation, SCA ’02, pages 167–174, 2002. ISBN 1581135734. → pages 4 D.T. Chen and D. Zeltzer. Pump it up: Computer animation of a biomechanically based model of muscle using the finite element method, volume 26. ACM, 1992. → pages 6 Nuttapong Chentanez and Matthias M´uller. Mass-conserving eulerian liquid simulation. In Proceedings of the 11th ACM SIGGRAPH/Eurographics conference on Computer Animation, pages 245–254. Eurographics Association, 2012. → pages 14 64 Michael Bradley Cline. Rigid body simulation with contact and constraints. PhD thesis, Citeseer, 2002. → pages 39, 40 S. L. Delp et al. OpenSim: open-source software to create and analyze dynamic simulations of movement. IEEE Trans Biomed Eng, 54:1940–1950, 2007. → pages 39 Franc¸ois Faure, J´er´emie Allard, Matthieu Nesme, et al. Eulerian contact for versatile collision processing. 2007. → pages 36 Nick Foster and Dimitri Metaxas. Realistic animation of liquids. Graph. Models Image Process., 58(5):471–483, 1996. ISSN 1077-3169. → pages 5 N. Galoppo, M.A. Otaduy, S. Tekin, M. Gross, and M.C. Lin. Soft articulated characters with fast contact handling. In Computer Graphics Forum, volume 26, pages 243–253. Wiley Online Library, 2007. → pages 52 Tolga G. Goktekin, Adam W. Bargteil, and James F. O’Brien. A method for animating viscoelastic fluids. ACM Trans. Graph., 23:463–468, August 2004. ISSN 0730-0301. doi:http://doi.acm.org/10.1145/1015706.1015746. URL http://doi.acm.org/10.1145/1015706.1015746. → pages 4 Patrick Hild. Numerical implementation of two nonconforming finite element methods for unilateral contact. Computer Methods in Applied Mechanics and Engineering, 184(1):99–123, 2000. → pages 8 K.R.S. Holzbaur, W.M. Murray, and S.L. Delp. A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control. Annals of biomedical engineering, 33(6):829–840, 2005. → pages 39 T.J.R. Hughes. The finite element method: linear static and dynamic finite element analysis. Dover Publications, 2000. → pages 5 John R. Humphrey, Daniel K. Price, Kyle E. Spagnoli, Aaron L. Paolini, and Eric J. Kelmelis. ¡title¿CULA: hybrid GPU accelerated linear algebra routines¡/title¿. 7705:770502–770502–7, April 2010. doi:10.1117/12.850538. URL http: 65 //proceedings.spiedigitallibrary.org/proceeding.aspx?doi=10.1117/12.850538. → pages 46 G. Irving, C. Schroeder, and R. Fedkiw. Volume conserving finite element simulations of deformable models. ACM Transactions on Graphics (TOG), 26 (3):13, 2007. → pages 5, 14 Kazufumi Ito and Karl Kunisch. Lagrange Multiplier Approach to Variational Problems and Applications. July 2008. URL http://dl.acm.org/citation.cfm?id=1413112. → pages 36, 44 K. Kamrin and J.C. Nave. An eulerian approach to the simulation of deformable solids: Application to finite-strain elasticity. arXiv preprint arXiv:0901.3799, 2009. → pages 5 K. Kamrin, C.H. Rycroft, and J.C. Nave. Reference map technique for finite-strain elasticity and fluid–solid interaction. Journal of the Mechanics and Physics of Solids, 2012. → pages 5 C. Lanczos. The Variational Principles of Mechanics. Dover Publications, 1986. → pages 10, 30, 35 Dongwoon Lee, Michael Glueck, Azam Khan, Eugene Fiume, and Ken Jackson. Modeling and simulation of skeletal muscle for computer graphics: A survey. Foundations and Trends in Computer Graphics and Vision, 7(4):229–276, 2012. → pages 6 E.H. Lee. Elastic-plastic deformation at finite strains. Technical report, DTIC Document, 1968. → pages 7 S.H. Lee, E. Sifakis, and D. Terzopoulos. Comprehensive biomechanical modeling and simulation of the upper body. ACM Transactions on Graphics (TOG), 28(4):99, 2009. → pages 2, 6 Sung-Hee Lee and Demetri Terzopoulos. Spline joints for multibody dynamics. In ACM Transactions on Graphics (TOG), volume 27, page 22. ACM, 2008. → pages 39 66 Michael Lentine, Mridul Aanjaneya, and Ronald Fedkiw. Mass and momentum conservation for fluid simulation. In Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’11, pages 91–100, New York, NY, USA, 2011. ACM. ISBN 978-1-4503-0923-3. doi:http://doi.acm.org/10.1145/2019406.2019419. URL http://doi.acm.org/10.1145/2019406.2019419. → pages 5, 31 David IW Levin, Benjamin Gilles, Burkhard M¨adler, and Dinesh K Pai. Extracting skeletal muscle fiber fields from noisy diffusion tensor data. Medical Image Analysis, 15(3):340–353, 2011a. → pages 7 D.I.W. Levin, J. Litven, G.L. Jones, S. Sueda, and D.K. Pai. Eulerian solid simulation with contact. ACM Transactions on Graphics (TOG), 30(4):36, 2011b. → pages 1, 5, 8, 10, 20, 36 Frank Losasso, Tamar Shinar, Andrew Selle, and Ronald Fedkiw. Multiple interacting liquids. ACM Trans. Graph., 25(3):812–819, July 2006. ISSN 0730-0301. → pages 4 H. Maas and T. G. Sandercock. Force transmission between synergistic skeletal muscles through connective tissue linkages. J. Biomed. Biotechnol., 2010: 575–672, 2010. → pages 39 GH Miller and P. Colella. A High-Order Eulerian Godunov Method for Elastic-Plastic Flow in Solids* 1. J. Comput. Phys., 167(1):131–176, 2001. ISSN 0021-9991. → pages 5 M.K. Misztal, K. Erleben, A. Bargteil, J. Fursund, B.B. Christensen, J.A. Bærentzen, and R. Bridson. Multiphase flow of immiscible fluids on unstructured moving meshes. In Eurographics/ACM SIGGRAPH Symposium on Computer Animation, pages 97–106. The Eurographics Association, 2012. → pages 5, 20 Richard M. Murray, S. Shankar Sastry, and Li Zexiang. A Mathematical Introduction to Robotic Manipulation. CRC Press, Inc., Boca Raton, FL, USA, 1st edition, 1994. ISBN 0849379814. → pages 30 67 Rahul Narain, Abhinav Golas, and Ming C. Lin. Free-flowing granular materials with two-way solid coupling. ACM Trans. Graph., 29:173:1–173:10, December 2010. ISSN 0730-0301. doi:http://doi.acm.org/10.1145/1882261.1866195. URL http://doi.acm.org/10.1145/1882261.1866195. → pages 6 P. Nardinocchi and L. Teresi. On the active response of soft living tissues. Journal of Elasticity, 88(1):27–39, 2007. → pages 7 V Ng-Thow-Hing, A Agur, and N McKee. A muscle model that captures external shape, internal fibre architecture, and permits simulation of active contraction with volume preservation. In INTERNATIONAL SYMPOSIUM ON COMPUTER METHODS IN BIOMECHANICS & BIOMEDICAL ENGINEERING (5.: 2001: Rome). Anais. Rome, 2001. → pages 13 Kiisa C Nishikawa, Jenna A Monroy, Theodore E Uyeno, Sang Hoon Yeo, Dinesh K Pai, and Stan L Lindstedt. Is titin a winding filament? a new twist on muscle contraction. Proceedings of the Royal Society B: Biological Sciences, 279(1730):981–990, 2012. → pages 7 S. Osher and R.P. Fedkiw. Level set methods and dynamic implicit surfaces. Springer-Verlag, 2002. ISBN 0387954821. → pages 5 T. Patterson, N. Mitchell, and E. Sifakis. Simulation of complex nonlinear elastic bodies using lattice deformers. ACM Transactions on Graphics (TOG), 31(6): 197, 2012. → pages 5 Alexei Y. Poludnenko and Alexei M. Khokhlov. Computation of fluid flows in non-inertial contracting, expanding, and rotating reference frames. Journal of Computational Physics, 220(2):678 – 711, 2007. → pages 6 Simon Premoˇze, Tolga Tasdizen, James Bigler, Aaron Lefohn, and Ross T. Whitaker. Particle-based simulation of fluids. Computer Graphics Forum, 22 (3):401–410, 2003. ISSN 1467-8659. doi:10.1111/1467-8659.00687. URL http://dx.doi.org/10.1111/1467-8659.00687. → pages 5 Karthik Raveendran, Chris Wojtan, and Greg Turk. Hybrid smoothed particle hydrodynamics. In Proceedings of the 2011 ACM SIGGRAPH/Eurographics 68 Symposium on Computer Animation, SCA ’11, pages 33–42, New York, NY, USA, 2011. ACM. ISBN 978-1-4503-0923-3. doi:http://doi.acm.org/10.1145/2019406.2019411. URL http://doi.acm.org/10.1145/2019406.2019411. → pages 6 Avi Robinson-Mosher, Tamar Shinar, Jon Gretarsson, Jonathan Su, and Ronald Fedkiw. Two-way coupling of fluids to rigid and deformable solids and shells. In ACM Transactions on Graphics (TOG), volume 27, page 46. ACM, 2008. → pages 37 Maurya Shah, Jonathan M. Cohen, Sanjit Patel, Penne Lee, and Fr´ed´eric Pighin. Extended galilean invariance for adaptive fluid simulation. In Proceedings of the 2004 ACM SIGGRAPH/Eurographics symposium on Computer animation, SCA ’04, pages 213–221. Eurographics Association, 2004. → pages 6 Juan C Simo and Robert L Taylor. Quasi-incompressible finite elasticity in principal stretches. continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering, 85(3):273–310, 1991. → pages 13 J. Stam. Stable fluids. In Proceedings of the 26th annual conference on Computer graphics and interactive techniques, pages 121–128. ACM Press/Addison-Wesley Publishing Co., 1999. → pages 5, 18 S. Sueda, A. Kaufman, and D.K. Pai. Musculotendon simulation for hand animation. ACM Transactions on Graphics (TOG), 27(3):83, 2008. → pages 2, 6, 40 Shinjiro Sueda, Garrett L. Jones, David I. W. Levin, and Dinesh K. Pai. Large-scale dynamic simulation of highly constrained strands. ACM Trans. Graph., 30:39:1–39:10, August 2011. ISSN 0730-0301. doi:http://doi.acm.org/10.1145/2010324.1964934. URL http://doi.acm.org/10.1145/2010324.1964934. → pages 5, 28 D. Sulsky, Z. Chen, and H.L. Schreyer. A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering, 118 (1-2):179–196, 1994. ISSN 0045-7825. → pages 5 69 J. Teran, S. Blemker, V. Hing, and R. Fedkiw. Finite volume methods for the simulation of skeletal muscle. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 68–74. Eurographics Association, 2003. → pages 6 J. Teran, E. Sifakis, S.S. Blemker, V. Ng-Thow-Hing, C. Lau, and R. Fedkiw. Creating and simulating skeletal muscle from the visible human data set. Visualization and Computer Graphics, IEEE Transactions on, 11(3):317–328, 2005. → pages 6, 13, 19, 22 Demetri Terzopoulos and Andrew Witkin. Physically based models with rigid and deformable components. Computer Graphics and Applications, IEEE, 8(6): 41–51, 1988. → pages 6 L. B. Tran and H. S. Udaykumar. A particle-level set-based sharp interface cartesian grid method for impact, penetration, and void collapse. J. Comput. Phys., 193(2):469–510, 2004. ISSN 0021-9991. → pages 5 J.A. Trangenstein. A second-order Godunov algorithm for two-dimensional solid mechanics. Computational Mechanics, 13(5):343–359, 1994. ISSN 0178-7675. → pages 5 Bin Wang, Franois Faure, and Dinesh K. Pai. Adaptive image-based intersection volume. ACM Trans. Graph. (Proc. SIGGRAPH), 31(4), 2012. → pages 20, 36, 53 Huamin Wang, Peter J. Mucha, and Greg Turk. Water drops on surfaces. ACM Trans. Graph., 24(3):921–929, July 2005. → pages 5 Q. Wei and D.K. Pai. Physically consistent registration of extraocular muscle models from MRI. In Engineering in Medicine and Biology Society, 2008. EMBS 2008. 30th Annual International Conference of the IEEE, pages 2237–2241. IEEE, 2008. → pages 7, 47, 55 Jeffrey A Weiss, Bradley N Maker, and Sanjay Govindjee. Finite element implementation of incompressible, transversely isotropic hyperelasticity. 70 Computer methods in applied mechanics and engineering, 135(1):107–128, 1996. → pages 13 Martin Wicke, Daniel Ritchie, Bryan M. Klingner, Sebastian Burke, Jonathan R. Shewchuk, and James F. O’Brien. Dynamic local remeshing for elastoplastic simulation. ACM Trans. Graph., 29:49:1–49:11, July 2010. ISSN 0730-0301. doi:http://doi.acm.org/10.1145/1778765.1778786. URL http://doi.acm.org/10.1145/1778765.1778786. → pages 6, 7, 23, 47 F.E. Zajac et al. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical reviews in biomedical engineering, 17(4):359, 1989. → pages 6 Yongning Zhu and Robert Bridson. Animating sand as a fluid. ACM Trans. Graph., 24(3):965–972, July 2005. → pages 5 71
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Eulerian finite volume method for musculoskeletal simulation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Eulerian finite volume method for musculoskeletal simulation and data-driven activation Fan, Ye 2013
pdf
Page Metadata
Item Metadata
Title | Eulerian finite volume method for musculoskeletal simulation and data-driven activation |
Creator |
Fan, Ye |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | This thesis describes a solid simulation method and its application to musculoskeletal simulation. The presented solid simulation method features Eulerian discretization and avoids mesh tangling during large deformation. Unlike existing Eulerian solid simulation methods, our method applies to elastoplastic material and volume-preserving material. To further increase the utility of Eulerian simulations for solids, we introduce Lagrangian modes to the simulation and present a new solver that handles close contact while simultaneously distributing motion between the Lagrangian and Eulerian modes. This Eulerian-on-Lagrangian method enables unbounded simulation domains and reduces the time step restrictions that often plague Eulerian simulation. We also introduce a framework for simulating the dynamics of musculoskeletal systems, with volumetric muscles and a novel muscle activation model. Muscles are simulated using the solid simulator developed and therefore enjoys volume preservation which is crucial for accurately capturing the dynamics of muscles and other biological tissues. Unlike previous work, in our system muscle deformation is tightly coupled to the dynamics of the skeletal system, and not added as an after effect. Our physiologically based muscle activation model utilizes knowledge of the active shapes of muscles, which can be manually drawn or easily obtained from medical imaging data. Finally we demonstrate results with models derived from MRI data and models designed for artistic effect. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2014-01-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0052197 |
URI | http://hdl.handle.net/2429/44645 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2013-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2013_fall_fan_ye.pdf [ 6.48MB ]
- Metadata
- JSON: 24-1.0052197.json
- JSON-LD: 24-1.0052197-ld.json
- RDF/XML (Pretty): 24-1.0052197-rdf.xml
- RDF/JSON: 24-1.0052197-rdf.json
- Turtle: 24-1.0052197-turtle.txt
- N-Triples: 24-1.0052197-rdf-ntriples.txt
- Original Record: 24-1.0052197-source.json
- Full Text
- 24-1.0052197-fulltext.txt
- Citation
- 24-1.0052197.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.24.1-0052197/manifest