University of Southern California Dissertations and Theses

Plant substructuring and real-time simulation using model reduction

(USC Thesis Other)

Plant substructuring and real-time simulation using model reduction

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Plant Substructuring and Real-time Simulation Using Model Reduction Yili Zhao Department of Computer Science Viterbi School of Engineering University of Southern California Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy August 2014 I would like to dedicate this thesis to my loving parents, Baochun Zhao and Tian Li. Acknowledgements Thank you, my advisor: Jernej Barbiˇ c. He is a genius, an expert in many fields. He gave me tremendous invaluable advices in academic research and life. He is a very gentle and nice mentor to work with, never ever yelling at his students (very important virtue for a PhD advisor!). It is my honor to have him be my PhD advisor. I will be nothing but grateful to him. Many thanks are given to my PhD thesis committee: Jernej Barbiˇ c, Gaurav S. Sukhatme, Ulrich Neumann, Igor Kukavica, and Stefan Schaal. I really appreciated all their construc- tive comments and feedbacks which helped me think critically and deeply into the research problems and greatly polish my thesis. Specially thanks Hao Li, a young, amazing, super talented, friendly, energetic, humorous professor. I received generous amount of help and advices from him! Sincerely thank you. Thank you Computer Science Department, USC. First of all, you are creative and innova- tive. You keep hiring so many top-notch scholars as faculties doing cutting edge research in broad areas of computer science, teaching us the knowledge, and mentoring us to the wonder- ful academic world. Second, you are friendly and supportive. Staffs including Lizsl De Leon Spedding, Michael Archuleta, Kusum Shori, Steve Schrader, Macy Yu really helped me a lot! Thank you. Last but not least, you are rich. As a PhD student here, I never worried about funding, just concentrated on study and research. Thank you USC, University of Southern California, (or rumor says University of Summer Construction) for your world class facilities and support. I love the Daland’s Swim Stadium (the main Aquatics venue at the 1984 Olympic Games). I swam a lot there, free style, 1 mile (nonstop) per day. Many thanks go to Computer Graphics Lab at Peking University. Thank you my teachers, iii Guoping Wang, Shihai Dong, Sheng Li. Without your immeasurable help, I could not fulfill my dream to pursue the PhD degree in a top university in the United States. Thanks to the fact that I do not have a girl friend for more than two decades since I entered the elementary school. I know it is a little bit pathetic, but it really saves me a lot of time, energy, and money, so that I can focus on study and doing research. Thanks my families: my parents, Baochun Zhao and Tian Li, my grandparents, aunts, uncles, and my cousins. No matter which path I choose, what adventure I take, you stand behind me. All of you are the source of my strength. I owe debts of gratitude to Rose Liu, Daniel Lin, David Lin, and Henry Lin, together with your parents. You are my families in Los Angeles. You made me feel I still stay at home, not in a foreign country. Thank you! Thank you, Lei Cheng, Chen Fan, Shengtao Zhou, Xin Huang. You are my brothers since childhood. I missed so much the old days when we had together. Thank you, Yao Sun, Dikang Gu, Bo Wang, Fenye Bao, Jian Zhao, Yi Sun, Ying Xiong, Jingjing Gu, my best friends during the college years. Thank you, Bin Lu, Chen Tang my best friends in Peking University. I missed the old time when we used to hang out for late-night snacks. I am much obliged to Wei Guan. You helped and guided me tremendously for those years in USC. You are like my elder-brother. Thank you, Xiaoyan Zhang. You are like my sister, always taking care of me. Thank you, Qiming Li. You are my little star. Thank you, Meihong Wang. I am always willing to listen to your opinions, exchange ideas with you. You are pretty smart and my best friend. Thank you, Jing Huang, the coach of ACM-ICPC team at USC, and my good roommate. We discussed a lot of algorithms in the spare time. Without these algorithms, I cannot publish my SIGGRAPH papers. Thank you, Yijing Li, Hongyi Xu, Guan Pang, Rongqi Qiu, Jiaping Zhao, Zhenzhen Gao, Lu Wang, Qianyi Zhou. I have a lot of good memories with all of you in the lab. Thank you, all my friends! The list could be incredibly long if I put all your names down. iv I cannot do that since it might be longer than my actual thesis. But please believe me you are always in my heart! Abstract This research is focusing on real-time, physically-based simulation of plants undergoing large deformations. To achieve this goal, we first propose a novel algorithm based on model reduc- tion and domain decomposition. It extends 3D nonlinear elasticity model reduction to open- loop multi-level reduced deformable structures. We decompose the input mesh into several domains, build a reduced deformable model for every domain, simulate each one separately, and connect domains using proper inertia coupling. This makes model reduction deformable simulations much more versatile: localized deformations can be supported without prohibitive computational costs, parts can be re-used and precomputation time can be shortened. Our method does not use constraints, and can handle large domain rigid body motion in addition to large deformations, due to our derivation of the gradient and Hessian of the rotation matrix in polar decomposition. We show real-time examples with multi-level domain hierarchies and thousands of reduced degrees of freedom. Then we design a pre-processor which takes a plant “polygon soup” triangle mesh as the only input and quickly pre-compute necessary data for the subsequent simulation. This tool breaks the ice for adoption of our multidomain dynamics simulator in practice. Our pre-processor is robust to non-manifold input geometry, gaps between branches or leaves, free-flying leaves not connected to any branch, small unimportant geometry (“debris”) left in the model, and plant self-collisions in the input configuration. Repeated copies (instances) of plant subparts such as leaves, petals or fruits can be automatically detected by our pre-processor. We enhanced our multidomain dynamics simulator to provide plant fracture, and inverse kine- matics to easily pose plants. It can simulate complex plants at interactive rates, subjected to user forces, gravity or randomized wind. We simulated over 100 plants from diverse cli- vi mates and geographic regions, including broadleaf (deciduous) trees and conifers, bushes and flowers. Our largest simulations involve anatomically realistic adult trees with hundreds of branches and over 100,000 leaves. Finally, we propose our future research in several directions including adding hierarchical in- stancing, collision detection and handling, etc. Contents Contents vii List of Figures x List of Tables xii Chapter 1 Introduction 1 1.1 Plant Simulation Using Multi-domain Dynamics . . . . . . . . . . . . . . . . 2 1.2 Plant Substructuring and Preprocessing . . . . . . . . . . . . . . . . . . . . . 2 1.3 Procedural Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Asynchronous Parallel Timestepping . . . . . . . . . . . . . . . . . . . . . . 3 1.5 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Chapter 2 Related Work 5 2.1 Plant Modeling and Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Finite Element Method and Model Reduction . . . . . . . . . . . . . . . . . 6 2.3 Domain Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Chapter 3 Multidomain Dynamics in Reduced Space 10 3.1 Background: Model Reduction and Domain Decomposition . . . . . . . . . . 10 3.2 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.5 Polar Decomposition Rotation Gradient . . . . . . . . . . . . . . . . . . . . 23 Contents viii 3.6 Appendix: System Force Integrals . . . . . . . . . . . . . . . . . . . . . . . 24 Chapter 4 Plant Preprocessing and Simulation 25 4.1 Plant Preprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2 Forests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.3 Interactive Plant Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.4 Fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.5 Motion of Leaves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 Chapter 5 Procedural Materials and Asynchronous Parallel Timestepping 51 5.1 Procedural Stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.2 Procedural Timestep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.3 Procedural Damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.4 Asynchronous Parallel Timestepping . . . . . . . . . . . . . . . . . . . . . . 58 5.4.1 Subdividing the Timestep . . . . . . . . . . . . . . . . . . . . . . . . 58 5.4.2 Parallel Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Chapter 6 Randomized Wind Force Field 63 6.1 Features of a Randomized Wind Force Field . . . . . . . . . . . . . . . . . . 63 6.2 Improved Perlin Noise Algorithm in 4D . . . . . . . . . . . . . . . . . . . . 64 6.3 Computing the Randomized Wind . . . . . . . . . . . . . . . . . . . . . . . 65 6.4 Applying Wind Forces to the Plant . . . . . . . . . . . . . . . . . . . . . . . 65 Chapter 7 System Overview 69 7.1 Botanical Preprocessor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 7.1.1 Command Line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 7.1.2 Preprocessor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 7.1.3 Useful Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 7.2 Botanical Simulator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 7.2.1 Simulation Control . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Contents ix 7.2.2 Frequency Tuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Chapter 8 Conclusion 81 Bibliography 84 List of Figures 1.1 Oregon white oak tree (Quercus Garryana) in strong randomized wind: . . . . 1 1.2 Anatomically realistic simulation of a peach tree (Prunus Persica) with fracture: 2 3.1 Domain decomposition: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2 Effect of system and interface forces: . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Fitting the best interface transformation. . . . . . . . . . . . . . . . . . . . . 14 3.4 Model reduction with a large number of localized degrees of freedom: . . . . 19 3.5 Instancing (Space Station): . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.6 Our method supports localized deformations. . . . . . . . . . . . . . . . . . 20 3.7 Similar trajectories, but different frequencies: . . . . . . . . . . . . . . . . . 21 3.8 C 0 embedding is effective with nearly-rigid interfaces: . . . . . . . . . . . . . 23 4.1 Preparing mechanical models ready for the simulation: . . . . . . . . . . . . 26 4.2 Branches (F), twigs (R1) and leaves (R2) . . . . . . . . . . . . . . . . . . . . 27 4.3 Instancing and anchors: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.4 Loop resolution: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.5 V oxel simulation meshes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.6 Hierarchy of a forest: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.7 Fir forest in randomized wind field: . . . . . . . . . . . . . . . . . . . . . . 39 4.8 Editing tree shapes using point constraints: . . . . . . . . . . . . . . . . . . . 40 4.9 Real-time Fracture: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.10 Leaf Motion: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.12 Representative subset of simulated models . . . . . . . . . . . . . . . . . . . 45 List of Figures xi 4.13 Space-time instantaneous (localized) external force followed by free vibration. 46 4.14 Palm tree (Cocos Nucifera) in strong randomized wind: . . . . . . . . . . . . 47 4.15 Gravity space: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.11 Visualization of domains in the oak tree: . . . . . . . . . . . . . . . . . . . . 48 4.16 Comparison to full simulation. . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.1 Avoidξ ’s minimum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.2 Asynchronous Timestepping: . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.3 Directed acyclic graph: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 6.1 Visualizing the Perlin wind field: . . . . . . . . . . . . . . . . . . . . . . . . 64 6.2 Applying Perlin wind to the plant: . . . . . . . . . . . . . . . . . . . . . . . 66 6.3 Grand fir (Abies Grandis) in Perlin wind: . . . . . . . . . . . . . . . . . . . 68 7.1 Plant preprocessor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 7.2 Select fronds and leaves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 7.4 Connecting the domains: . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 7.5 Resolving Loops: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 7.3 Select root domain: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 7.6 Botanical Simulator: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 7.7 Simulation control panel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 7.8 Frequency tuning panel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 List of Tables 3.1 Simulation statistics (multidomain dynamics) . . . . . . . . . . . . . . . . . 19 4.1 Timings for each pre-processing step . . . . . . . . . . . . . . . . . . . . . . 30 5.1 Simulation statistics (asynchronous parallel timestepping) . . . . . . . . . . . 62 Chapter 1 Introduction A large fraction of our world is covered by vegetation. Botanical environments are both di- verse and very common; therefore, they are crucial for special effects, games and virtual reality applications. The goal of this thesis is to design a real-time, physically-based simulator to com- pute large-deformation dynamics of plants. Figure 1.1: Oregon white oak tree (Quercus Garryana) in strong randomized wind: the first fully mechanically simulated, botanically realistic and accurate, adult tree, anywhere in sci- ence, 871 branches, 120,000 leaves, 2,360,868 triangles, 9,520 reduced DOFs, simulation fps: 6 Hz 1.1 Plant Simulation Using Multi-domain Dynamics 2 Figure 1.2: Anatomically realistic simulation of a peach tree (Prunus Persica) with fracture: Peaches fall from the tree swaying in the randomized wind. 299,707 triangles, 237 branches, 3,556 twigs, 18,536 leaves, 330 fruits, 2,950 reduced DOFs, 7 hierarchy levels, simulation fps: 20 Hz 1.1 Plant Simulation Using Multi-domain Dynamics We first propose a novel algorithm based on domain decomposition and model reduction which is capable of simulating complex, botanically accurate models with large deformation at inter- active rates. Given a volumetric mesh of a plant, we decompose the mesh into several pieces called domains, build a reduced deformable model for each domain, simulate them separately, and connect them using inertia coupling. This makes model reduction deformable simulations much more versatile: localized deformations can be supported without prohibitive computa- tional costs, parts e.g. twigs and leaves can be re-used and pre-computation time is shortened. Our method does not use constraints, and can handle large domain rigid body motion in ad- dition to large deformations, due to our derivation of the gradient and Hessian of the rotation matrix in polar decomposition. We show real-time examples with multi-level domain hierar- chies and thousands of reduced degrees of freedom (figure 1.1 and figure 1.2). 1.2 Plant Substructuring and Preprocessing Existing domain decomposition methods do not address the problem of how to robustly and quickly pre-process many plants in the presence of imperfections in input geometry, which is a key challenge for the adoption of such methods in practice. We design a plant pre-processor 1.3 Procedural Materials 3 which is robust to non-manifold input geometry, gaps between branches or leaves, free-flying leaves not connected to any branch, small unimportant geometry (“debris”) left in the model, and plant self-collisions in the input configuration. We remove loops in the domain graph using a new user-assisted algorithm to select a minimum spanning tree in a general undirected graph. Our novel domain graph creation algorithm, instancing, and spanning tree selection procedure apply to any domain decomposition plant simulation method, including those that do not employ model reduction [96]. Our system supports plants represented using triangle meshes and alpha-masked billboards and exploits multi-core computation. Real-time fracture is supported, enabling our plants to shed leaves or drop fruits (Figure 1.2). We also present an approach to perform physically-based inverse kinematics, enabling the user to adjust the geometry of existing plants by dragging plant vertices. We simulated over 100 plants from diverse climates and geographic regions, including many broadleaf (deciduous) trees and conifers, bushes and flowers. Our system can simulate both simple plants and plants as complex as entire anatomically realistic adult trees with several hundreds of branches and over 100,000 leaves (see Figure 1.1). 1.3 Procedural Materials Tuning simulation parameters is quite critical in order to generate vivid and plausible results. For geometrically complex plants, such tuning process is quite challenging since it is too com- plicated to do it manually in practice. In this thesis, we propose a series of procedural methods to automatically tune simulation parameters including material stiffness, timestep, and damp- ing coefficients for botanical objects with large structural and geometrical complexity. These tuning techniques are generic and can be applied to any modal-based FEM simulation. 1.4 Asynchronous Parallel Timestepping We have designed an asynchronous parallel timestepping algorithm to timestep the domains more efficiently. Since natural frequency is different for each domain, it inspires us to integrate 1.5 Thesis Overview 4 each domain at its own pace as opposed to using a fixed timestep for the entire object. This avoids the shortcomings from a synchronous integrator, in which a small, locked time step is carefully chosen in order to prevent from numerical instabilities due to the high stiffness variations encountered at different domains. Our total running time to produce animations of the same length (same total physical time) is therefore shorter than in the synchronous method. By exploiting the fact that no loop exists in the tree topology, we extend simulation to multi- cores, introduce buffers to decrease serializing dependencies, and optimize the structure of our multi-domain dynamics algorithm. We demonstrate that such parallelization greatly enhances simulator’s performance (5x). 1.5 Thesis Overview The rest part of this thesis is organized as follows. In Chapter 2, we provide a comprehensive survey of related work. This survey includes 3 parts: plant modeling and simulation, finite element method (FEM) and model reduction, and domain decomposition (substructuring). In Chapter 3, we present a novel method which makes model reduction adaptive in space. We call it multi-domain dynamics algorithm. In this algorithm, we first decompose the deformable object into several components (called domains), then pre-process the reduced dynamics of each domain separately, and finally couple the domains using inertia forces. In Chapter 4, we demonstrate a robust plant preprocessor, which takes a static plant triangle mesh as the only input and pre-computes the mechanical models ready for the simulation use. In Chapter 5, we propose a series of procedural methods to tune simulation parameters for complex botanical objects. These properties include material stiffness, timestep, damping coefficients, etc. We design an asynchronous parallel timestepping method based on the multi-domain dynamics. In Chapter 6, we describe how we generate a randomized wind field in our simulator. We adopt the improved Perlin noise [77] which is widely used in texture synthesis and create a 4D randomized wind field. In Chapter 7, we provide a brief overview of our botanical simulation pipeline which includes two main parts: botanical preprocessor and simulator. In Chapter 8, we conclude our research and list several directions for future work. Chapter 2 Related Work 2.1 Plant Modeling and Simulation Plant modeling has a long history in computer graphics and we refer the reader to good survey work [23]. There are many strategies to create and edit plant geometry [78, 79]. Lindenmayer systems (L-systems) are probably the most well-known system for plant modeling and simu- lation [57, 67, 79]. Image-based techniques produce 3d plant models from photographs [80], and can be augmented with sketch-based interfaces to guide the plant development [36]. Plant geometry can also be guided by specific environment conditions [78]. While the early plant modeling methods were purely procedural [14, 72], user-assisted plant modeling [58] has been demonstrated to scale and produce realistic models of complex plants, as seen in the commer- cial Xfrog system [101]. We use Xfrog models extensively in our work, as well as models from other sources, such as the SpeedTree [37] system, and other online 3d model stores [95]. Our method works with “triangle soup” input meshes and is agnostic of the specific modeling approach. Perhaps the simplest approach to animating plants is to not use physics at all; but drive the deformations kinematically using a stochastic wind [100]. Animations can also be created us- ing pre-recorded motion graphs [42, 105], optionally combined with stochastic motion [104]. Physical simulation, however, has the advantage that it can provide dynamics, large motion, secondary motion and easier runtime control. A common physical approach is to model 2.2 Finite Element Method and Model Reduction 6 branches as rigid rods, and connect them with angular springs [83, 98, 103]. Such an approach is usually augmented with a proper randomized wind model [73], and can be accelerated using level of detail techniques [10]. In nature, branches are flexible and their bending has been modeled using one-dimensional flexible beams [32, 34], and using the three-dimensional solid Finite Element Method (FEM) [59, 96]. Three-dimensional FEM simulations have the advan- tage that they can work with “polygonal soup” input geometry, simply by converting them to volumetric meshes. They can easily support volume preservation and spatially-varying ma- terial properties. They also automatically incorporate branch thickness (thicker branches are harder to bend) and non-straight (crooked) undeformed branch geometry. Our simulator builds upon the 3D solid Finite Element Method (FEM). It augments it with frame-aware domain de- composition to support large deformations, modularity, interactive design and fracture, and employs model reduction to increase the speed of computation. 2.2 Finite Element Method and Model Reduction Deformable object simulation is a well-studied problem in computer graphics. We review approaches for interactive FEM and model reduction; please see [68] for a general survey. FEM simulations with complex geometry do not run at interactive rates. Interactivity can be achieved using multi-resolution geometric constructions [17, 22, 31], employing co-rotational elasticity [19, 65], the multigrid method [27], or by coarsening of meshes and their mate- rial properties [46, 69]. In our work, we employ model reduction (see [87] for a survey), and demonstrate that whenever the object can be decomposed into natural components, this can provide deformation-rich real-time simulations. The common theme in model reduction is to project the original state-space onto a low-dimensional subspace to arrive at a (much) smaller system having properties similar to the original system [51]. The idea has been em- ployed in several disciplines outside of computer graphics, such as control theory [28], elec- trical circuit simulation [55], computational electromagnetics [88], microelectromechanical systems [38], fluid simulation [33, 50, 60] and nonlinear elasticity [2, 51, 70]. Many pow- erful reduction techniques have been designed, both for linear (in many cases time-invariant) 2.2 Finite Element Method and Model Reduction 7 systems [5], such as balanced truncation [64] and Krylov subspaces [30], and nonlinear sys- tems [53, 81]. Despite this progress, model reduction of large-scale, high-dimensional (thou- sands of degrees of freedom) nonlinear systems is still a relatively new field [11]. In com- puter graphics, model reduction of nonlinear systems has been used for fast simulation of deformable solids [3, 8, 44, 49, 62] and fluids [94], and for fast control of such systems [7]. For deformable FEM offline simulations, Kim and James [47] showed how to apply online model reduction to adaptively replace expensive full simulation steps with reduced steps. One drawback of these systems has been that the reduction basis is global in space. Wicke and colleagues [99] extended Treuille’s fluid reduction method to several inter-connected reduced domains. Their approach is similar in spirit to ours, but works for fluids and does not directly apply to nonlinear elasticity. Our open-loop domain structures are related to recursive algo- rithms for articulated rigid body structures [26]. Most of these algorithms were developed for rigid objects. Motivated by the applications in robotics, Featherstone’s algorithm has been extended to deformable objects [15, 18, 86], to small-deformation simulations [18, 86], and deformable rods [12]. Our algorithm, however, applies to 3D solid deformable objects with irregular interfaces undergoing large deformations. Application of model reduction to linear systems are common [76], and offer several advan- tages such as rapid simulation rates, as well as easy runtime control over the simulated and ren- dered level of detail. In the context of plant simulation, Stam [91] used a modal basis to com- pute noise in the frequency domain, and Diener [24] used a wind projection basis to increase computation speed. Low-dimensional (three of less) linear modal simulations on individual, fully decoupled, branches modeled as Euler-Bernoulli beams were presented in [32, 34], as well as methods to tune the models to match recorded real tree motion. These approaches produce quality high-frequency motion of trees, e.g., leaves and branches rapidly fluttering in the wind, and are as such complementary to our nonlinear model reduction method. Because they employ linear models, they cannot, however, support large deformation plant dynamics, characteristic of non-wooden (herbaceous) plants or trees blowing in strong winds. Research on how to quickly simulate plants undergoing large deformations, especially the ones with large structural and geometrical complexity, appears very limited. In our research, we develop 2.3 Domain Decomposition 8 a robust system which allows us to quickly pre-process many complex botanically accurate models and launch them in our interactive simulation program. 2.3 Domain Decomposition The idea to decompose a deformable object into several interconnected components (domains), each of which can be simulated separately, is well-known in science. It is usually referred to as domain decomposition [93], or, especially in case of repetitive geometry, also substructur- ing [25, 71]. Many such methods do not employ reduction, but merely divide the object so that each domain can be assigned to a different processor node, or data for repetitive substructures can be reused [6, 82]. Perhaps the simplest form to add reduction to domain decomposition is to compute the static equilibrium of each domain, under x,y,z perturbations of each inter- face vertex, and then restrict the domain deformations to a linear combination of those shapes (static condensation [92]). With complex geometry, however, interfaces themselves can be high-dimensional, leading to a large number of basis vectors and slow simulation times. Al- ternatively, one can compute the linear vibration modes of each structure, under the boundary condition that the interfaces are held fixed [21]. This method, called component mode synthe- sis [84], has been popular with simulations of deformations of superstructures, most notably in aerospace engineering (e.g., airplanes, space satellites) [75]. These previous methods, how- ever, only simulated small, linearized deformations of each domain. It is not straightforward to extend them to large deformations because the resulting large interface rotations seemingly require modes to rotate, which invalidates precomputation. In our paper, we show how these obstacles can be avoided, yielding a component mode synthesis method supporting (1) large deformations within each domain, and (2) large (finite) rotations of the domain interfaces. This makes substructuring much more useful in computer graphics applications requiring large de- formations. Domain decomposition for deformable models has also been previously applied in computer graphics, but only for small domain deformations and with running times dependent on the number of domain and interface vertices. For example, a linear quasi-static applica- tion using Green’s functions has been presented in [41], whereas Huang and colleagues [35] 2.3 Domain Decomposition 9 exploited redundancy in stiffness matrix inverses to combine linear FEM with domain decom- position. Because plants typically consist of well-defined subparts with only limited interaction, domain decomposition is a natural approach to tackling plant motion complexity. Twigg and Kaˇ ci´ c- Alesi´ c [96] simulated deformable objects consisting of many parts, including trees, by gluing them together using the Procrustes transform. Their method simulates each deformable part in the full space. Such simulations can be accelerated by orders of magnitude using model re- duction, which makes it possible to timestep complex plants at interactive rates. Recently, the combination of domain decomposition and model reduction has received significant attention in the literature [9, 48, 102]. Our research [9], for example, employed gradients of rotation matrices computed in polar decomposition, for structures that do not have loops in the domain graph. Kim and James [48] efficiently simulated deformations of characters using inter-domain spring forces. Recent work of [102] uses modes obtained from unit displacements of interface vertices, and inertia modes, combined with modal warping. Our method does not require a skeleton, pre-existing motion, or well-defined domain inter- faces and it supports fracture, inverse kinematics, frequency tuning, and bottom-up meshing. Different from [48], a global volumetric mesh is not necessarily needed, as each domain is meshed individually, enabling component re-use and modularity. Chapter 3 Multidomain Dynamics in Reduced Space 3.1 Background: Model Reduction and Domain Decompo- sition Model reduction is a popular method for deformable model simulation, mainly because it can approximate complex physical systems at a low computational cost. The key idea of model reduction is to project the high-dimensional equations of motions to a suitably chosen low-dimensional space where the dynamics have properties similar to the original system, but can be timestepped much more quickly [51]. Real-time projection-based model reduction for deformable objects has, however, suffered from an important limitation: the reduction basis is global in space and time. Such bases require a large number of modal vectors to capture local deformations. More importantly, because nonlinear modal elasticity requires implicit integration for stability, and because all global basis vectors overlap in space, each timestep requires (at least) solving a ˆ r× ˆ r dense linear system costing O(ˆ r 3 ), where ˆ r is the number of basis vectors. In practice, this has limited real-time nonlinear reduced simulations to less than (approximately) one hundred degrees of freedom [3]. In this Chapter, we present an approach to make model reduction adaptive in space, by decomposing the deformable object into several components (the domains, see Figure 3.1). We pre-process the reduced dynamics of each domain separately, and then couple the domains 3.2 Kinematics 11 using inertia forces. Assuming a decomposition free of loops, the resulting system supports large deformation dynamics both globally and locally within each domain. For the geomet- rically nonlinear FEM material model, the resulting nonlinear system can be timestepped at rates independent of the underlying geometric or material complexity. With exact reduced internal force evaluations on d domains with r degrees of freedom each, the running time of one timestep of our method is O(dr 4 )≪ O(ˆ r 4 ), for ˆ r= dr, and could be further decreased to O(dr 3 ) using approximate reduced forces [3]. The idea of decomposing a deformable object for efficient simulation has been previously extensively explored in the engineering community, usually under the names of domain de- composition and substructuring. However, previous methods either did not pursue reduction in each domain, or limited the domains to small deformations. Our method is related to the well-known Featherstone’s algorithm for linked rigid body systems, but differs from it by sim- ulating large deformations involving large interface rotations, combined with model reduction. Like the Featherstone’s algorithm supports kinematic chains of arbitrary length, we make no assumptions on the depth of hierarchies. We approximate subtree inertia using mass lumping, which gives us fast and stable real-time large deformations rich in local detail. Our method supports instancing: a single object can be pre-processed once and replicated many times sav- ing on precomputation and runtime costs. Rigid domains (r= 0) are supported, and can be arbitrarily mixed with reduced-deformable domains. The method also supports unanchored objects undergoing free-flight motion. 3.2 Kinematics Our method uses reduction to simulate geometrically nonlinear FEM deformations of a 3D vol- umetric mesh. We show examples with linear tetrahedral elements and trilinear voxel (cube) elements. Let a volumetric mesh be decomposed into d connected and mutually disjoint sets of elementsD 1 ,D 2 ,...,D d (the domains, see Figure 3.1). The common surface where two domains i and j meet is called the interface, I i j . Although in principle an object could be de- composed arbitrarily, domain decomposition is most effective when the domains form a natural 3.2 Kinematics 12 root (b) (a) (d) (c) root Figure 3.1: Domain decomposition: (a) domains, (b) domain graph, (c) 48 domains (space sta- tion), (d) domain graph (space station). decomposition of the object, such as, for example, separating the space station modules, panels and antennas into separate domains (Figure 3.1). In such cases the interfaces are often small or deform mostly rigidly (see Figure 3.2, right), which we exploit to define our low-dimensional kinematic model. We first form the domain graph, where each domain is one node, and nodes are connected if they share a common interface (see Figure 3.1). We assume that there are no cycles in the domain graph (graph is a tree) but place no restriction on maximum node degree. We direct the graph by picking one domain as the root node, which then uniquely directs all edges by traversing the domains from the root to the leaves. We set the root domain to the do- main rooted to the ground for fixed objects, or a central/largest domain for free-flying objects, which tends to minimize tree depth. Figure 3.2: Effect of system and interface forces: Left: Pulling on the small stem causes large secondary motion (bloom, leaves). Right: relative deviation of the interface transformation A i j from rotation R i j (Frobenius norm, for all the 11 interfaces). The deviations are small. We simulate each domain as a nonlinear reduced-deformable object, with its own specific 3.2 Kinematics 13 reduced basis U i of size r i ≥ 0. All our reduced models are precomputed under the boundary condition that the domain is held fixed at the interface to the parent, which is consistent with the interface rigidity assumption. The deformation of each domain in its local frame of reference, is given as a linear combination of the modes of domain i, u i = U i q i , where u i contains the 3D deformations of all the vertices of domain i, for the reduced-coordinate vector q i ∈R r i . Deformations u i are expressed in a local frame of referenceF i of each domain, which we define below. Frames are necessary because parts of the mesh may undergo large rotations, whereas modal models are poorly suited to represent large rotations. At first, it may seem possible to avoid frames by including all affine transformations into each basis U i . However, the remaining vectors in U i (the non-rigid deformations) would then need to be rotated synchronously with the domain rotation at runtime, leading to a time-dependent basis and a significant additional computational cost. To define the frames, we first collect all individual vectors q i into a global vector q∈R r 1 +...+r d . The frame computation then proceeds from the tree root to the leaves. FrameF 0 is the world coordinate frame for fixed objects, and the global rigid body motion frame for free-flying objects. For each child domain j of domain i, we define frame F j as the best fitting frame to the interface I i j . We do so by specifying its position x i j ∈R 3 and rotation R i j ∈R 3× 3 , relative to frame F i , and expressed in the coordinate axes of frame F i , by fitting the best rigid transformation that transforms vertices of interface I i j from their rest positions to current positions given by q i (see Figure 3.3). Let v 1 i j ,...v m i j be the vertices of domain i that are on the interface to child domain j. We can weight the vertices according to the surface area (or mesh volume) locally belonging to each vertex, arriving at weights w 1 i j ,...w m i j . In domain i, each vertex k deforms according to a 3× r i submatrix of U i , denoted by U k i j . We make x i j track the centroid of I i j : x i j = 1 ∑ m k=1 w k i j m ∑ k=1 w k i j (X k i j +U k i j q i )≡ ˆ x i j + ˆ a i j q i , (3.1) where X k i j is the rest position (inF i ) of vertex v k i j , and the matrix ˆ a i j ∈R 3× r i can be precom- puted. In order to fit the rotation, we need to align the interface vertices to their deformed 3.2 Kinematics 14 positions as best as possible using a rigid transformation (after subtracting centers). We do so by first computing the covariance matrix [66] A i j (q i )= B i j (q i )C − 1 i j ˆ R i j ≡ ˆ R i j + r i ∑ ℓ=1 A ℓ i j q ℓ i , where (3.2) B i j (q i )= m ∑ k=1 w k i j (X k i j − ˆ x i j )+(U k i j − ˆ a i j )q i X k i j − ˆ x i j T , (3.3) C i j = m ∑ k=1 w k i j X k i j − ˆ x i j X k i j − ˆ x i j T , (3.4) and the fixed matrices A ℓ i j ∈R 3× 3 can be precomputed. Here, ˆ R i j is the rotation of interface I i j in the rest configuration, relative to rest frame F i . We then perform polar decomposition to extract the best-fitting rotation R i j : A i j (q i )= R i j (q i )S i j (q i ), (3.5) REST DEFORMED D j D i Figure 3.3: Fitting the best interface transformation. where S i j is a symmetric 3× 3 matrix. We note that rigid transformations cannot be made to fit deformable interfaces in general, leading to discrepancies in the domain meshes at the interface. In practice, however, the interface transforma- tion are often very close to rigid and the discrepancies are small (Figures 3.2). If interface vertices lie in the same plane in the rest configuration, matrix C i j becomes degenerate. Let 0≤ λ 1 ≤ λ 2 ≤ λ 3 be the eigenvalues of C i j . Ifλ 1 <ελ 3 (we use ε = 0.001), we modify the precomputation of A ℓ i j by adding an extra point at ˆ x i j +ηN, where N is the eigenvector forλ 1 (inter- face normal). The value η = p ελ 3 − λ 1 increases the smallest eigenvalue of C i j toελ 3 , and gave stable planar interfaces in our examples. Once the root frameF 0 and transformations(x i j ,R i j ) are known for all interfaces I i j , we can easily compute world-coordinate expressions for all framesF i , i= 0,...,d− 1. We note that the framesF i are completely determined by q andF 0 ; i.e., the frames are not separate inde- 3.2 Kinematics 15 pendent simulation parameters. An alternative would be to keep them separate and simulate the combined frame-deformable system, but doing so would require constraints to keep the domains connected, turning the system into a differential-algebraic equation and leading to standard problems of constraint drift. Our construction, in turn, avoids constraints. In the next section, we shall derive the dynamics that govern q, using a mass-lumped formulation running in time linear in the number of domains. This is similar to Featherstone’s algorithm, but for flexible objects undergoing large deformations. For dynamics, it is necessary to compute linear velocity v i j , linear acceleration a i j , angular velocity ω i j , and angular acceleration α i j of each frameF j , relative toF i and expressed in F i , based on q i , ˙ q i , ¨ q i . Linear quantities are simply v i j = ˆ a i j ˙ q i and a i j = ˆ a i j ¨ q i . To compute the angular quantities, first differentiate Equation 3.2: ˙ A i j = r i ∑ ℓ=1 A ℓ i j ˙ q ℓ i , ¨ A i j = r i ∑ ℓ=1 A ℓ i j ¨ q ℓ i . (3.6) Next, we compute the first and second derivatives of the polar decomposition rotation matrix R i j with respect to q i . We derived a general formula for ˙ R(t) and ¨ R(t), where R(t) is the rotation in polar decomposition of a 3× 3 matrix A(t)= R(t)S(t) that depends on a scalar parameter t∈R (usually time, but can be any parameter): G= tr(S)I− S R T ∈R 3× 3 , ω = G − 1 2skew(R T ˙ A) ∈R 3 , (3.7) ˙ R= ˜ ωR, ˙ S= R T ( ˙ A− ˙ RS), ¨ R= ˜ ˙ ωR+ ˜ ω 2 R, (3.8) ˙ ω = G − 1 2skew R T ( ¨ A− ˜ ω ˙ A) − tr( ˙ S)I− ˙ S R T ω . (3.9) Here, ˜ ω denotes the 3× 3 skew-symmetric matrix corresponding to a vector ω ∈R 3 , i.e., ˜ ωx=ω× x for all x∈R 3 . Similarly, skew(A) denotes the unique skew-vectorω∈R 3 so that ˜ ω =(A− A T )/2. The derivation of Equations 3.7 - 3.9 is given in Section 3.5. Evaluation of ω, ˙ R, ˙ S, ˙ ω, ¨ R requires solving the 3× 3 nonsymmetric linear system given by matrix G. Because this system is nonsingular whenever A is nonsingular, which is the case for interfaces that deform mostly rigidly, computing G − 1 is stable. We can now apply Equations 3.7-3.9 3.3 Dynamics 16 to A i j (t)= A i j (q i (t)) and its time derivatives as given by Equations 3.2 and 3.6, yielding ω i j andα i j . We note that such matrix decomposition gradients have been previously explored for singular value decomposition [61, 96]. With SVD, singularities in the decomposition gradient occur whenever two singular values are equal, e.g., even if A is identity, and would, unlike polar decomposition, require additional treatment. 3.3 Dynamics We now give the equations of motion for q, under the kinematic model of Section 3.2. These equations simulate the coupled motion of all domains. Each domain follows the equation M i ¨ q i + D i ˙ q i + f int i (q i )= f ext i + f sys i + ∑ j is child of i f itf i j , (3.10) where M i is the reduced mass matrix (constant matrix), D i = D i (q i ) is the reduced damping matrix, and f int i (q i ) are the reduced nonlinear internal elastic forces of domain i. The terms f sys i , f itf i j , f ext i contain the reduced system forces due to motion of ancestor domains, interface forces of child domains of i arising due to subtree inertia, and external forces, respectively. Equation 3.10 is standard in model reduction; it is obtained by projecting a full (geometrically nonlinear) FEM deformable model to a chosen basis U i . We use the modal derivative basis [8] because the computation is automatic and does not require any presimulation. We address non-inertiality in by computing non-inertial frame of reference system forces f sys i due to the motion of frame i as imposed by its parent domain. System Forces: Equation 3.10 is expressed in frame F i which is non-inertial (accelerates through time). An observer rigidly attached toF i can correctly simulate deformations if she adds the resulting system forces f sys i to the equations of motion of her domain (see, e.g. [40]). Let X be a material point in frameF i . The world-coordinate velocity and acceleration of X, 3.3 Dynamics 17 expressed in the coordinate frameF i , equal [85] v(X)= v i +ω i × X+ ˙ X (3.11) a(X)= a i +α i × X+ω i × (ω i × X)+ 2ω i × ˙ X+ ¨ X, (3.12) where v i ,ω i ,a i ,α i are the world-coordinate velocity, angular velocity, acceleration and angular acceleration of frameF i , respectively, expressed in the frameF i . The system forces are dis- tributed volumetrically throughout the domain. When projected to the low-dimensional space of each domain, they are f sys i =− Z D i ρ i (X)U T i (X)a(X)dV, (3.13) where ρ i (X) is mass density at X, U i ∈R 3× r i are the spatially-varying modes, and a(X) is acceleration at X. Theℓ-th component of f sys i ∈R r i , forℓ= 1,...,r i , can be expanded to f sys iℓ = W 1T iℓ a i − W 2T iℓ α i +(W 3 iℓ + q T i W 4 iℓ )||ω i || 2 − − (ω T i W 5 iℓ + 2 ˙ q T i W 6 iℓ )ω i − q T i W 6 iℓ α i − (ω i ω T i ) : r i ∑ p=1 W 7p iℓ q p i , (3.14) for constant precomputable coefficients W j iℓ (Section 3.6). Notation A : B denotes component- wise matrix dot product. The evaluation of f sys i requires O(r 2 i ) flops, and is fast in practice. Interface Forces: We model reduced interface forces as f itf i j =− M i j ¨ q i + f 0 i j , for (3.15) M i j = m j ˆ a T i j ˆ a i j , f 0 i j = ˆ a T i j R i j f ext j − (3.16) − m j a i +α i × x i j +ω i × (ω i × x i j )+ 2ω i × v i j , where m j and f ext j are the total mass and net sum of external forces (expressed inF j ) in the subtree rooted at j (call itD j ). Note that m j is constant and precomputable, but could vary 3.3 Dynamics 18 at runtime, e.g., if domains fracture, or are replaced (modularity). Equation 3.15 is similar to Featherstone’s recursive terms. The terms M i j and f 0 i j model the mass inertia of D j , by assuming that the mass is lumped at I i j . The term f 0 i j also models the effect (to domain i) of external forces applied inD j . This approximation gives stable motion and keeps the system matrix symmetric, and we found it reasonable for examples with limited domain graph depth. Integration: Equation 3.15 transforms Equation 3.10 into M i + ∑ j M i j ¨ q i + D i ˙ q i + f int i (q i )= f ext i + f sys i + ∑ j f 0 i j . (3.17) We timestep Equation 3.17 in time linear in the number of domains. The algorithm first con- structs the frames for all domains, then timesteps each domain using semi-implicit Newmark integration for stability [8], and finally updates relative frame kinematics (see Algorithm 1). The frame of the root domain can be the world-coordinate frame if the root domain is fixed to the ground, or its frame can be floating and affected by the forces and torques of the subdo- mains for free-flying objects. Algorithm 1: Multidomain reduced dynamics 1 Procedure Simulation timestep Input: values of ˆ q (k) =(q, ˙ q, ¨ q) and framesF (k) at timestep k, external forces f ext,(k+1) at timestep k+ 1, timestep size∆t Output: values of ˆ q (k+1) andF (k+1) at timestep k+ 1 2 begin 3 for all domains i, leaves to root do Assemble f ext i ,m i 4 for all domains i, root to leaves do 5 Computeω i ,a i ,α i ,F (k+1) i (Eq. 3.11, 3.12) 6 Compute f sys i (Eq. 3.14) 7 for all child domains j of i do 8 Compute M i j , f 0 i j (Eq. 3.16) 9 Do a reduced dynamics timestep (Eq. 3.17), yielding ˆ q (k+1) i 10 for all child domains j of i do 11 Use ˆ q (k+1) i to compute x i j ,R i j ,v i j ,ω i j ,a i j ,α i j 3.4 Results 19 3.4 Results Example vol-vtx vol-el rend-vtx rend-tri d ˆ r D pre rt:core rt:StVK rt:total S u i = U i q i fps flower (tet) 2,713 7,602 6,675 12,868 12 240 3 0.6 min 22 % 78 % 6.5 msec 1 0.43 msec 69 Hz SIGGRAPH (tet) 18,945 83,753 10,463 20,934 15 160 8 4.3 min 21 % 79 % 1.7 msec 1 2.5 msec 83 Hz dragon (tet) 46,736 160,553 38,625 77,250 40 454 5 5.1 min 45 % 55 % 1.6 msec 3 7.0 msec 40 Hz space station (voxel) 219,058 107,556 177,691 248,521 48 921 4 0.8 min 25 % 75 % 4.5 msec 3 25.5 msec 14 Hz oak tree (tet) 262,363 626,734 578,801 838,704 1435 11,972 5 1.0 min 52 % 48 % 29 msec 1 12.2 msec 5 Hz Table 3.1: Simulation statistics for #volumetric and rendering mesh vertices, elements and tri- angles (vol-vtx, vol-el, rend-vtx, rend-tri), #domains (d), total # of reduced DOFs (ˆ r), tree depth (D), precomputation time (pre), #simulation steps per graphics frame (S), constructing deforma- tions for rendering (u i = U i q i , once per frame), frame rate (fps), simulation step time (rt:total), and its breakdown in terms of timestepping reduced dynamics of individual domains (rt:StVK) and the rest (rt:core, including polar decomposition, gradient and Hessian, system and interface forces, frames). Machine specs: Intel Core i7-980X, 6-Core, 3.33 GHz, 10 GB memory. Only a single core was used. Figure 3.4: Model reduction with a large number of localized degrees of freedom: Left: non- linear reduced simulation of an oak tree, all leaves are deformable domains (41 branches (r= 20), 1394 leaves (r= 8), d= 1435 domains, ˆ r= 11,972 total DOFs) running at 5 fps. Right: simulation detail. In our first example we show a complex oak tree (41 branches, 1394 leaves) deformed in the wind and undergoing large deformations (see Figure 3.4). This example greatly uses sub- structuring, as the leaves all use a single mesh (rotated and translated into the proper place). The leaf mesh is pre-processed and reduced only once. The copies use pointers to the single datastructure, leading to short precomputation times (Table 3.1). Similarly, the 41 branches only belong to 5 distinct classes, translated, rotated and scaled (uniformly) into their place. 3.4 Results 20 Figure 3.5: Instancing (Space Station): The repeated panel copies are instanced, and can all bend independently, as can all major structure parts. Scalings can be handled efficiently by observing that for a uniform scaling factor s, the ba- sis matrix does not change, whereas the frequency spectrum scales by 1/s. We also pursued instancing in the space station example (Figure 3.5). Figure 3.6: Our method supports localized deformations. End-effector domains often deform most due to largest system forces. The number of modes in each domain can be chosen automatically, by setting a total num- ber of modes ˆ r, and assigning to each domain the number of modes proportional to its number of elements. One good choice is a logarithmic distribution: r i ∝ log(#elements(i)) (dragon ex- ample, Figure 3.6). The SIGGRAPH example (Figure 3.6) demonstrates that our method sup- 3.4 Results 21 ports free-flying motion. This is achieved by pre-processing a “free-fly” reduced deformable model [8] for the root domain, and then using the external forces to integrate the rigid body motion for the entire object. Our method supports rigid domains (r i = 0), so it can combine rigid and deformable objects: the 8 letters of “SIGGRAPH” are connected by 7 rigid links. Figure 3.7: Similar trajectories, but different frequencies: Top and middle: Bird’s-eye view on the trajectories of two flower vertices: top of primary (A) and secondary bloom (B). Same scale used for A and B. Rest and extreme poses are indicated. Ground truth (G)=solid black, interface lumping (L1)=dashed green, center of mass lumping (L2)=dashed red. Bottom: x-dof of A versus time. Accuracy: Figure 3.7 compares our method to an unreduced single-domain geometrically nonlinear FEM simulation (ground truth). We provide a comparison to two variants of our method: (L1) interface lumping (Equation 3.16), and (L2) center-of-mass lumping where the mass of the subtree is lumped not at the interface, but at the center of mass of the subtree. Method L2 can be implemented by adding additional terms to Equation 3.16, guided by Equa- tion 3.12. The ground truth (G) was computed offline and is 55x slower than L1 and L2. All three simulations use the same material parameters, timestep, fixed vertices, and the initial condition: a velocity aligned with the first eigenmode of the entire mesh, sufficient to bend the flower stem by about 45 degrees. We found that the three methods give similar trajectories, but differ in oscillation frequencies: L1 and L2 gave 2x and 1.5x higher lowest natural frequency than G, respectively. In general, reduced simulations of solids lack detailed DOFs, resulting 3.4 Results 22 in small increases in natural frequency (“artificial stiffening”). With multidomain reduced dy- namics, however, the increase is largely due to mass lumping, similar to how a pendulum with mass lumped close to the pivot oscillates at a higher frequency than if the mass was distributed further away. L2 matches G more closely than L1 because lumping at the subtree center better approximates the actual mass distribution. There is a pyramid of methods that model the subtree inertia progressively better, at the cost of additional implementation complexity. All these methods take the form of Equation 3.15, but differ in how M i j and f 0 i j are computed. All examples in Table 3.1 use method L1, which we found to be the simplest approach producing stable, deformation-rich results. In L1, each do- main feels the total weight of the attached subtree, but not the rotational or deformable inertia. Variant L2 improves the accuracy for physically long domains. In some examples, however, we observed that L2 suffers from instabilities; we attribute this to the quickly time-varying matrices M i j in method L2. At a significant additional implementation complexity, one could compute correct, non-lumped subtree inertia for kinematic chains with arbitrary tree depth (cf. [85]). Such M i j and f 0 i j would fully parallel Featherstone’s algorithm, but may require modifications to the integrator to maintain stability. With uniform material parameters, small mesh parts (e.g., protrusions) vibrate at higher fre- quencies than the rest of the mesh, which manifests as little or no deformations in those regions even with the ground truth (unless poked explicitly). The frequency content of each domain can be linearly scaled at runtime without redoing the precomputation, simply by scaling the domain’s precomputed reduced internal forces. In the flower example (Figure 3.2), we ad- justed the frequencies of the leaves so that the leaves gave interesting large deformations when pulling on the stem. Similarly, we made the horns, tail, spike and mouth of the dragon softer than the rest of the mesh, to cause larger deformations in those regions (Figure 3.6). Rendering: We render triangle meshes embedded into volumetric meshes. Although the volumetric mesh for the entire object is a manifold mesh in the rest configuration, the domains slightly separate at the interfaces under deformation. At first, we anticipated this to be prob- lematic – however, we found that it can be very easily handled with techniques similar to those 3.5 Polar Decomposition Rotation Gradient 23 employed in discontinuous Galerkin FEM [45]. We can establish C 0 continuity simply by av- eraging the two copies of each interface vertex (see Figure 3.8), whereas C 1 continuity could be achieved via moving least square (MLS) embeddings [45]. Once a consistent volumetric mesh deformation is computed, it is transferred to the embedded triangle mesh using barycen- tric interpolation. Because our bases are local in space, computing the vertex deformations via equation u i = U i q i only involves matrices U i with a small number of columns, leading to small memory footprints and fast computation. Figure 3.8: C 0 embedding is effective with nearly-rigid interfaces: (a),(c): individual domain meshes, for two representative deformations. Domain gaps (black) are generally small; most are sub-pixel size. (b),(d): Mesh deformations with C 0 continuity. 3.5 Polar Decomposition Rotation Gradient Let A = A(t) be a 3× 3 matrix that depends on a scalar parameter t ∈R. For any t, one can perform polar decomposition, A(t)= R(t)S(t), where R is orthogonal and S is symmetric positive semi-definite. Fix t 0 ∈R. We first shift the problem by defining B(t)= R T (t 0 )A(t)= R T (t 0 )R(t) S(t). (3.18) Because R T (t 0 )R(t) is identity for t = t 0 , its derivative at t = t 0 must be a skew-symmetric matrix ˜ ω for some ω ∈R 3 . By differentiating Equation 3.18 by t, and setting t = t 0 , one obtains ˙ B(t 0 )= R T (t 0 ) ˙ A(t 0 )= ˜ ωS(t 0 )+ ˙ S(t 0 ). (3.19) 3.6 Appendix: System Force Integrals 24 We now replace t 0 with t, and apply the skew operator (Section 3.2) to both sides of Equa- tion 3.19. This causes the symmetric term ˙ S to drop and yields 3 linear equations for the 3 components ofω (Equation 3.7). In order to derive ˙ ω and ¨ R, differentiate both sides of Equa- tion 3.7 with respect to t (note that G and ω depend on t). After rearranging, one obtains a 3× 3 linear system for ˙ ω (Equation 3.9). 3.6 Appendix: System Force Integrals The constants of Equation 3.14 are integrals over each domainD i : W 1 iℓ =− Z D i ρU ℓ i dV∈R 3 , W 2 iℓ = Z D i ρ ˜ XU ℓ i dV∈R 3 , (3.20) W 3 iℓ = Z D i ρU ℓT i XdV∈R, W 4 iℓ = Z D i ρU T i U ℓ i dV∈R r i , (3.21) W 5 iℓ = Z D i ρU ℓ i X T dV∈R 3× 3 , W 6 iℓ = Z D i ρU T i f U ℓ i dV∈R r i × 3 , (3.22) W 7p iℓ = Z D i ρU ℓ i U pT i dV∈R 3× 3 , (3.23) where U ℓ i ∈R 3 is theℓ-th column of mode U i (X)∈R 3× r i . Chapter 4 Plant Preprocessing and Simulation Previous computer graphics research on plants has focused on plant geometry creation and appearance (rendering), as well as efficient simulation. Simulation of complex plants is chal- lenging, however, because plant meshes are typically designed for rendering, not simulation. In this chapter, we demonstrate how to robustly and quickly pre-process complex plants in the presence of imperfections in the input geometry, for subsequent fast physically based simula- tion. Because plants naturally decompose into their constituent parts (branches, twigs, leaves, etc.), this makes botanical objects fit for our multi-domain dynamics simulator described in the previous chapter. Such authoring of simulation-ready plants augments and completes sim- ulation, by making it easy to apply our model-reduction based simulator to general, complex plant models. 4.1 Plant Preprocessing Given an input triangle mesh of a plant, we present an efficient and easy to use user-assisted pipeline to organize the plant into domains, create the domain hierarchy, create the simulation meshes, and pre-process reduced models (Figure 4.1). The end result is a plant mechanical model that can be rapidly simulated using our model reduction simulator (previous chapter). All of the steps of our pipeline except the specific model reduction pre-processing in Step 10 apply generally to plant domain decomposition, and could be used with other domain decom- 4.1 Plant Preprocessing 26 Figure 4.1: Preparing mechanical models ready for the simulation: (a) Triangle mesh of an European larch (Larix Decidua). This mesh was designed for rendering and is not suitable for simulation. (b) Simulation meshes: the tree is divided into different domains (16,549 domains, shown in different colors). We compute a mechanical model for each domain seperately. The total number of DOFs is 9,825. (c) Real-time simulation (using wind), using the preprocessed mechanical models. Simulation is running at 30 frames per second. Rendering is at 5 frames per second. 4.1 Plant Preprocessing 27 position methods [48, 96]. Our procedure starts with a triangle mesh of a plant. We make no assumptions on mesh topology or connectivity and support arbitrary “triangle soups” which may contain cracks, T-vertices or duplicated triangles. Such plant models are very common in practice. The various plant parts need not be properly connected to each other in terms of sharing vertices; for example, it is sufficient if adjacent branches simply collide with each other slightly at their common intersection (see Figure 4.2, right). Our system supports “bill- boards”, i.e., texture-mapped (usually simple) triangle meshes with transparency, commonly used to model leaves, fronds and smaller branches (twigs). We now explain our pre-processing pipeline. Steps performed by the user are marked as (U), whereas fully automated steps are marked as (A). Table 4.1 analyzes the time needed for each of the steps. Figure 4.2: Branches (F), twigs (R1) and leaves (R2) Branches can simply “sink” into each other in the input geometry (right). Step 1: (U) Organize input mesh into domains: Given a “polygon soup” mesh, the poly- gons must be grouped into domains. Each of our domains is characterized by the user as one of three types: (i) F, (ii) R1 or (iii) R2 (see Figure 4.2). Domains F are flexible, and typically in- corporate meshed, non-billboard 3D geometry, e.g., the trunk, branches, or flowers. Domains of type Rx are rigid, and are in practice often instanced. They are typically used for fruits, billboard twigs, fronds, leaves, small decorative geometry or even unwanted geometry left in the model by artists (“debris”; e.g., small unconnected triangles). We use two levels R1 and R2 because in some plants, the artist intended billboard domains to be parented to other billboard domains, e.g., conifer billboard needles attached to billboard twigs. If a single level R1 was 4.1 Plant Preprocessing 28 employed, some domains may be parented to flexible domains that are too far away, which can cause neighboring domains to separate at runtime. A typical example of the decomposition is to assign all branches into F, twigs into R1 and leaves into R2 (Figure 4.2). The user is free to deviate from such guidelines, however. For example, with some models, leaves are modeled as detailed triangle meshes, in which case they may be considered of type F. Every triangle must be assigned into exactly one domain, and each domain is assigned one of the F,R1,R2 types. We employ a user interface similar to that in, say, Maya, where the user can select triangles or domains, and show / hide / add / subtract / delete / merge them. Domains can be selected with the mouse and then tagged as either F, R1 or R2. In some models, polygons are pre-grouped by the artists into individual logical parts, e.g., each leaf is a separate domain already in the input, in which case the domains must only be selected, and their type identified. Often, however, the parts of the same kind are grouped together, e.g., all leaves are initially in one domain, and must therefore be separated into individual leaf domains. Therefore, another operation that we support is to break an existing domain into connected components, where two triangles are considered connected if they share a vertex. We compute the connected components using the union-find datastructure [20]. When breaking domains into components, the user can choose to impose a rule that all triangles in the domain may use at most one texture image; otherwise, the domain is broken further, into one domain for each texture image. Such a rule is useful with instancing, but is rarely needed because most models already satisfy this requirement as is. Figure 4.3: Instancing and anchors: (a) One-to-one texture map with transparency. (b) User- selected anchor points (in purple). 4.1 Plant Preprocessing 29 Step 2: (U) Instancing: Many plants consist of repeated parts. The instances are translated, rotated and sometimes scaled copies of each other, e.g., replicated leaves or flower petals. We automatically identify instances as follows. We assume that the triangle mesh of each instanceable domain is texture-mapped with a single image, with a one-to-one texture map (Figure 4.3, a). In particular, vertices cannot have different(u,v) values if they appear in more than one triangle, and triangles cannot “fold over” or cross each other in the(u,v) space. In our model databases, we did not encounter texture maps that would not be one-to-one. We identify instances as follows. We first inspect the texture image and the number of vertices; domains that do not match in both, are not instanced copies. For each vertex i of the first domain, we then find the vertex j=P(i) of the second domain whose texture coordinates are closest to those of vertex i, by performing a nearest neighbor search in the (u,v) space. If the mapP is not one-to-one and onto (a permutation), or the distance to the nearest neighbor is greater thanε = 10 − 3 for any vertex, we deem the domains not instanced. Finally, we check that the triangle mesh topology is same for both domains, i.e., if vertices i, j,k form a triangle in the first domain, so must verticesP(i),P( j),P(k) in the second domain, and vice-versa. If topology is the same, domains are deemed instanced copies of each other. At first, we attempted to use instancing where we did not seek the permutationP, but simply checked the (u,v) distance between vertex i of first domain and vertex i of second domain. This approach did not work well because modeling packages sometimes arbitrarily re-order triangle mesh vertices before exporting them. Note that our definition of instancing relies on the (u,v) space and therefore permits arbitrary translation, rotation and scaling of the object in world-space. All the domains that are instanced copies of one another are placed into an instance set. A single domain is chosen as a representative instance; we choose the one that appears first in the input mesh. For each instance set, we also ask the user to select an “anchor point” on the representative instance (see Figure 4.3, b). The user selects the anchor point in world coordinates by clicking with the mouse on the model, upon which we determine and store the corresponding(u,v) coordinate. The anchor point will be used to determine the instance transformation, and also in Step 7 to assign a parent to each Rx domain. The anchor point should typically be selected at the end of the botanical piece, e.g., on leaf’s stem (petiole) next to the attachment point to branch (see 4.1 Plant Preprocessing 30 Example 1 2 3 4 5 6 7 8 10 total dahlia 16 40 0.2 0 3 27 0.01 0.45 141 228 jasmine 29 14 0.6 20 3 63 0.09 1.13 124 255 white pine 36 3 7 10 3 60 3.50 17 170 311 broadleaf 7 10 1.4 0 3 269 1.14 3.3 136 431 Table 4.1: Timings for each pre-processing step, in seconds, for four representative models (shown in Figure 4.12), including user and computer time. Four users were involved in the experi- ments. Timings for steps 9 and 11 are not shown to save space, but are included in totals; they are less than 0.7 sec in all examples. Figure 4.3, b). Sometimes, billboard domains are not designed to attach to anything, but are free-floating in great numbers, e.g., clusters of needles on a conifer tree (Figure 4.3, b1). In such cases, the user should select the center of the billboard polygon. Next, for each instance from an instance set, we determine the linear transformation that best aligns the vertices of the representative instance to this instance, using shape matching [66]. Optionally, the user can force the linear transformation into a rotation, computed using polar decomposition. It is often preferable to use linear transformations, however, because the leaves are often scaled in size to add variety. The instance transformation is stored to disk, and used at runtime to properly transform the domain. Only a single mesh must be stored for each instance set. Step 3: (A) Computing the F domain graph: Our system then automatically builds the domain graph for the F domains. The nodes of the graph are the F domains, and two nodes are connected if the two triangle meshes intersect in the undeformed configuration. We determine the graph edges using collision detection, accelerated by bounding volume hierarchies and spatial hash tables [56]. Because we only need to perform collision detection once (on static shapes), the domain graph construction only takes a few seconds at most, even for our most complex examples. Initially, we attempted to use collision detection on volumetric meshes (computed using voxelization) to create our hierarchy. Although such an approach avoided gaps between domains, it created many spurious graph edges, which greatly complicated span- ning tree selection in Step 6. 4.1 Plant Preprocessing 31 Step 4: (U) Connecting the F domain graph: The graph computed in step 3 may not be connected. In practice, we encountered such situations in about 25% of all models. Most often, this occurs for one of two reasons: (1) There is a small gap between two domains, often visually (nearly) invisible and unimportant for rendering, or (2) the domains are “debris”: small pieces inadvertently left in the model by the artist, often invisible inside branches. We let the user connect the graph as follows. We compute the connected components (using union- find). The user can then scroll among the components, select arbitrary two domains in arbitrary components, and connect them. At any time of this process, she can recompute the graph and the connected components. In most cases, when the graph was initially not connected, the total number of components was less than 10. We encountered a few cases of “debris” with approximately 50 components, which we deleted one by one. When visualizing the connected components, we sort them base on their size. At any moment, the user can opt to simply delete the remaining components, if their size is deemed insignificant. Step 5: (U) User selects the root domain: Plants in nature are rooted. Although most models come pre-oriented so that the Y-axis is up, this is not guaranteed, so the user has to specify the root domain. Typically, this is the stem or the trunk. We provide the domain with the largest diameter as the initial suggestion. Step 6: (U) User-assisted resolution of loops: Although in principle the output of Step 4 should have no cycles (a tree), this is not always the case in practice. For example, an artist may have accidentally left a branch colliding with another branch, forming a cycle (see Figure 4.4, top). There are three very common causes of cycles: collisions of non-neighboring branches in the input (very common with complex trees), “Y-bifurcations” where a branch splits into two branches and all three meshes collide with each other, forming a 3-cycle, and petals on flowers (Figure 4.4, D). We found cycles to be common in commercial plant model libraries. They must be removed so that a tree hierarchy can be computed. Because the computer lacks con- text, it is difficult to remove cycles automatically, e.g., the small transverse branch in Figure 4.4 (magnified in A), may as well be connected to either of the main branches. A human, how- ever, can look at the branch and recognize its intended direction. Therefore, we designed an 4.1 Plant Preprocessing 32 Figure 4.4: Loop resolution: A: a 7-cycle. B: an example of a 3-cycle where a computer may make a mistake and that requires user intervention. Because both smaller branches collide with the bigger branch, computer initially suggests an incorrect loop-breaking graph edge (center), whereas the correct edge is shown on the right. C: the minimum spanning tree selection algorithm. (1): Initial minimum spanning tree (green) and redundant edges (red). The edge to be resolved next is marked by “add”. (2): The loop. User decides to cut the edge marked with “X”. (3): The minimum spanning tree after deleting the edge. (4): The hierarchy after the remaining two redundant edges were processed. D: Examples of loops in input geometry: “Y”-bifurcation and flower petals. 4.1 Plant Preprocessing 33 algorithm for user-assisted removal of loops. Our algorithm applies generally to any problem where a quality minimum spanning tree must be selected out of the many minimum spanning trees of an undirected graph with cycles. We first compute an initial candidate minimum span- ning tree, using a breadth-first traversal starting from the root domain. If the graph has no cycles, we are done; otherwise, the edges which are not in the minimum spanning tree are added to the set of redundant edges, prioritized by their breadth-first order (Figure 4.4, C1). We then let the user resolve these edges, one by one. For each redundant edge, we add it to the minimum spanning tree, and therefore exactly one cycle appears in the resulting subgraph S (Figure 4.4, C2). The cycle is computed by traversing the minimum spanning tree from each of the two vertices of the redundant edge towards the root of the tree, until a common ancestor is detected, thereby detecting a cycle. The cycle is then visualized to the user, by coloring the cycle domains in a golden color (as in Figure 4.4, top). The user is then asked to break the cycle by removing exactly one edge. The user does so by scrolling through the cycle (as in Figure 4.4, top; the two node domains that are the endpoints of the edge are shown in red and pink), and selecting the edge to be removed. After the edge removal, the working graph S is a tree again, and the redundant edge set has shrunk by one (Figure 4.4, C3). This process is repeated until the redundant set becomes empty, i.e., the working graphS is a tree and all the redundant edges were resolved. We then orientS to form a tree hierarchy, starting from the root and proceeding to the leaves (Figure 4.4, C4). Because redundant edges are priori- tized by their breadth-first traversal order, the user resolves cycles closer to the root first. In our tree database, virtually all plants initially had a non-tree domain graph, and required user intervention. Most loops, however, occur due to “Y” bifurcations, and can be resolved very quickly. For small / moderate examples, these were often the only loops. The vast majority of redundant edge sets that we have encountered in our examples had less than 100 edges. A few large examples, such as large trees, had approximately 500 redundant edges. Even for the most complex trees, the redundant edge set removal was manageable and was completed within minutes of user time. We were able to significantly shorten the user time by implement- ing an auto-focus feature where the camera automatically focuses on each loop as the user is scrolling through the loops. 4.1 Plant Preprocessing 34 Step 7: (A) Add domains Rx to the domain tree: We now assign parents to domains Rx. Domains R1 are always parented to F domains, whereas domains R2 can be parented by R1 or F. We use anchor points to determine the parent. The anchor position is computed using the inverse of the texture map, by finding the world-coordinate location on the domain mesh whose texture coordinates are (u,v), where (u,v) are the texture coordinates selected on the representative instance by the user in Step 2. Such anchor position computation was very robust in practice. We then perform a nearest neighbor search, seeking the nearest triangle in all domains in F (and also in R1 for R2) to the anchor position; the closest domain becomes the parent. Such an assignment is robust to cracks and can accommodate (intentionally) “floating” domains that are common with billboarding. For example, leaves on trees in many models are simply accumulated in close proximity to give a space-filling perception, and are far from any branch. The nearest neighbor query is accelerated using a bounding volume hierarchy (we use axis-aligned bounding boxes), and only takes a few seconds, even for complex models (Table 4.1). Figure 4.5: Voxel simulation meshes. Left: root domain. Middle, Right: two representative domains. Note that the meshes for the different domains are completely independent and do not need to meet in a common interface. Adjacent petals can thus be meshed independently without any difficulty, even though they collide in the input triangle mesh (middle, right). Fixed vertices are shown in red, and need not be vertices of the parent volumetric mesh; some are even outside the parent mesh. Step 8: (A) Build simulation meshes: We build the volumetric mesh for each domain by voxelizing the domain’s triangle mesh [39] (Figure 4.5). We chose this approach because it is completely automatic and supports arbitrary (potentially ill-formed) input “polygon soup” 4.1 Plant Preprocessing 35 geometry. Alternatively, one could employ automatic tetrahedral meshes such as [52]. The user specifies the maximum voxelization resolution (it defaults to 100 in our system) for the expanded bounding box of the input triangle mesh (we use expansion factor of 1.2× ). The voxelization resolution of each domain is then set automatically, by assigning to each domain a resolution proportional to the longest edge of its bounding box. Using such an adaptive reso- lution, the meshes of different domains grade approximately uniformly across the entire plant (see Figure 4.5). We note that the voxel meshes for the individual domains are separate meshes; they need not topologically connect to meshes of other domains using any well-defined inter- face. Such meshing flexibility is possible in domain decomposition simulators that do not need the connectivity requirement, such as [9, 96]. Our “bottom-up” mesh generation approach has the advantage that an entire simulation mesh never needs to be constructed. It is also very modular. Input triangle mesh collisions of neighboring pieces (e.g., petals on a flower) are not a concern, because each piece is simply meshed separately (see Figure 4.5). If a new piece needs to be added, it can be pre-processed in isolation and added to the domain graph. In [9], domain decompositions were created by forming a global tetrahedral mesh and then manually subdividing it into domains. Although such a ”top-down” approach results in well-defined in- terfaces, it requires careful manual work to produce quality interfaces. For example, the choice of whether to place a tet into the left or the right domain at the interface can affect the bend- ability of short branches. Also, top-down approaches require special handling in the presence of collision in the input triangle mesh; otherwise, the colliding branches or petals will simply be welded in the global volumetric mesh. Step 9: (A) Assign fixed vertices: Fixed vertices serve as connection points of a domain to its parent. The root domain is rooted to the ground, either by fixing all vertices that are below a user-provided height, or by manually selecting its fixed vertices. Fixed vertices for the remaining domains are set automatically as follows. For each parent-child pair in the hierarchy, we constrain the child vertices that are located close to the parent. Specifically, we detect child voxels that intersect parent voxels, and constrain all the vertices in those child voxels. Because we always constrain all eight vertices of a voxel, we avoid degenerate planar sets of fixed 4.1 Plant Preprocessing 36 vertices. The intersection test is accelerated using a bounding volume hierarchy. Such an assignment handles both the case where the child voxels are smaller than parent voxels and may be completely inside a parent voxel, and the case where child voxels are large and may subsume parent voxels. It can happen that all the vertices of a domain are deemed fixed; e.g., with small, in-grown, branches close to a bigger branch; in such cases, we declare the domain to be rigid. When no fixed vertices are detected, the domain is also made rigid. Such cases are rare in our data, but sometimes occur with domains not properly connected to the rest of the plant, e.g., small branches added for decoration in the tree crown. Note that for our domain decomposition approach, it is not necessary to establish a well-defined interface between the two domains; the child fixed vertices may even be outside of the parent mesh (Figure 4.5, middle, right). Also note that in botanical systems, the interfaces between branches and/or leaves are often small in surface area, and their bending can be neglected. Namely, the state- of-the-art papers in plant simulation go even further and assume that all the parts (the branches) are completely decoupled [32, 34]. In our work, the domains are coupled, with the assumption that the interface deformation is small. Step 10: (A) Compute low-dimensional simulation basis and pre-process reduced dynam- ics: Given the fixed vertices, we compute linear modes, their large-deformation correction (modal derivatives [8]), and the simulation basis. By default, we use the first 10 linear modes, and a 20-dimensional basis. We use a default (and tunable), spatially uniform, mass density ρ= 1000kg/m 3 , Young’s modulus E= 10 6 N/m 2 and Poisson’s ratio ofν= 0.45, correspond- ing to a fairly incompressible material. We then precompute a geometrically nonlinear FEM reduced model for each domain [8]. Each of these models is a compact, low-dimensional representation of the FEM dynamics of each domain. It supports large deformations, and can be timestepped rapidly, in microseconds. The domains are coupled to each other as de- scribed in Section 3.2. The entire process is automatic. In order to save space and increase speed, we use an adaptive number of modes per domain, by pre-processing smaller bases for smaller domains. The user sets the maximum size of a domain simulation basis r max (we use r max = 20). For a domain with e elements, the number of modes is then computed as 4.2 Forests 37 r = max(5,⌊r max log(e)/log(e max )⌋), where e max is the maximum number of elements in a domain (see also [9]). Step 11: (A) Tune Young’s modulus: After the reduced models are pre-processed, they have correct stiffness levels relative to each other. However, the global stiffness scale is arbitrary; one can globally scale ρ and E with an arbitrary constant without affecting the basis, and the precomputed reduced polynomials only scale by a constant [8]. Because we are using a FEM method, it would be possible to set the material properties for wood, stems, leaves, etc. Although progress on measuring elastic material properties has been made [13], the parameters for plants are typically not available. Therefore, we opt to assign default material values, pre-process the reduced models, then exploit the fact that natural vibration frequencies scale linearly under a linear global scale of Young’s modulus. We scale Young’s modulus so that the lowest natural vibration frequency of the root domain becomes 1 Hz, i.e., E ′ = E/ f 2 0 , where E is the default Young’s modulus, and f 0 is the lowest natural frequency of vibration of the root domain, determined using a sparse eigenvalue solver [54]. The user can then further scale the spectrum as needed, to make the plant stiffer or softer. Tuning is very fast at runtime; it takes less than one second even for complex models. 4.2 Forests It is easy to expand our simulator to simulate a forest (Figure 4.7). The key idea is to form a virtual root to accommodate the forest. It becomes the root of our hierarchy (Figure 4.6). This virtual root is simply a box, large enough so all the botanical objects can be rooted. It is rigid and fixed in the world coordinate frame, thus, no simulation is needed. All the botanical objects on the top are pre-processed separately. Replicated copies of the same plant form an instancing group and only the representative of this group needs to be pre-processed. We continue to spell-check how to handle transformation of an instanced copy. Scaling If the mesh is scaled simultaneously in x, y, z by factorα, the reduced deformation basis must be scaled byα − 3/2 . Note that, in this case, the tangent stiffness matrix K is scaled 4.3 Interactive Plant Design 38 Figure 4.6: Hierarchy of a forest: (a) (b) Hierarchies for two plants. (c) A virtual root (shown in red) is created to accommodate the forest. It becomes the root of the hierarchy. Replicated plants (plant 1 shown in yellow) form an instancing group. by α, the mass matrix is scaled by α 3 , and the frequency is scaled by1/α. These scalings happen automatically in our system. User may have to tune the stiffness properties in the simulator to meet the application needs. Rotation The reduced deformation basis must be rotated accordingly. Often with plants, however, rotation is not needed. Translation Translate the plant to the correct position. No special operation is needed for the pre-processed data. 4.3 Interactive Plant Design Because our method is fast and modular, it is possible to use it as an interactive shape editor for plants. The user can easily delete or add new parts at runtime, either by duplicating existing plant parts, or importing parts from a pre-processed library of plant parts. The editing process is local and interactive. The mass, stiffness, position and orientation of each branch, twig or leaf, as well as a linear scale (geometric size), can all be adjusted interactively at runtime, with immediate physically based simulation feedback to the user, without any additional pre- 4.3 Interactive Plant Design 39 Figure 4.7: Fir forest in randomized wind field: It consists of 20 firs, each of which has 81,228 triangles, 9,077 domains, 4,858 DOFs. Entire scene is simulated at 2 frames per second, including u= Uq construction of vertex deformation. Rendering is at 1 frames per second. processing. It is possible to randomize these choices, creating an arbitrary number of variations of the same plant. In our interactive editor, we can select an arbitrary domain, and linearly scale the stiffness (Young’s modulus) of all of its elements by any factorα> 0. It can be shown that under such a scaling, the reduction basis does not change [8], whereas the reduced forces scale by α; therefore, the scaling is instant, and there is no need to maintain an explicit volumetric mesh or its element material properties. Because stiffness is proportional to the square of the lowest natural vibration frequency, such scaling makes it possible to cause a branch to oscillate faster or slower. For cinematic effect (to show more secondary motion), we sometimes found it useful to make smaller branches oscillate more slowly than dictated directly by physics. We achieve this by scaling every domain with a factor α = β d , where β > 0 is a constant and d is the depth of domain in the hierarchy. Because our method is fast, the user can tune β interactively with immediate feedback; typically, values close toβ = 0.1 were producing good results. In addition, in our system the user can also manipulate the plant using inverse kinematics-like 4.3 Interactive Plant Design 40 Figure 4.8: Editing tree shapes using point constraints: The tree shape was adjusted to avoid collision with the house. The tree manipulated points are shown in red. Bottom row shows the tree without the leaves. 45 branches, 125 twig billboards, 1154 leaf billboards, 8 msec of simulation for one graphical frame. 4.3 Interactive Plant Design 41 handles. Using such a tool, it is possible to constrain and drag plant vertices, while the rest of the plant automatically re-adjusts using physics to a good-looking, minimal strain energy configuration. For example, a tree can be made to lean in a certain direction or avoid external objects (Figure 4.8). The user can also use it to resolve any unwanted collisions in the rest configuration, or simply re-adjust the plant shape to increase scene geometric variety. In our IK tool, the user can select or deselect an arbitrary number of vertices (IK handles). As the user drags the mouse, a three-dimensional force is applied to the active IK vertex, whereas the remaining IK vertices are kept fixed to their positions using linear springs. The mouse force is applied in the “screen plane”, i.e., plane orthogonal to the view direction and cutting through the manipulated vertex. The force magnitude is proportional to the number of pixels traveled by the mouse since the beginning of the drag. We use our standard real-time dynamics solver for such IK manipulation; no special code is required. Because the model normally undergoes dynamic motion, we must employ some mechanism to quickly stabilize the motion to a limit equilibrium configuration. One approach is to perform a Newton-Raphson iteration to seek the model equilibrium under the IK linear spring forces (a static solver [7, 63]). However, such a solver often suffers from a high linear system condition number, which has lead to instabilities in our experiments. We found much better results by simply using a high level of stiffness-proportional Rayleigh damping, in a dynamic simulation. Although the model does not reach the equilibrium configuration instantly, with proper gains and damping the convergence is rapid and stable. A single stiffness gain for all the linear springs was sufficient in our examples. Because we already injected sufficient damping into our simulator, linear springs did not need any additional damping. We note that we first attempted to perform inverse kinematics by enforcing exact user control over the handle positions. This was performed by solving, at every timestep, an optimization problem that minimized the total strain energy of the plant subject to the exact IK constraints. This approach did not work very well in practice. Because the handles must be specified in the three-dimensional space which can be difficult to visualize on a 2D screen, it was easy for the user to accidentally command unreasonable handle positions, inconsistent with plant dimensions. As the solver was trying to meet the constraint exactly, the plant would stretch 4.4 Fracture 42 unnaturally, which led to vibrations and instabilities. Instead, when positions are enforced via springs, the solver has the ability to selectively “yield” on each constraint as needed. This tends to produce much smoother, natural-looking shapes. We were able to tune the IK stiffness gains to reach both good output shapes and minimal deviation of the constrained vertices from their prescribed positions (Figure 4.8). After the user is satisfied with the shape, she can save it to disk, and re-process the reduced models with respect to the new rest shape. For small/moderate edits where the old reduced basis is still sufficient, a re-process is not necessary; in this case, one can simply compute the reduced forces for the new shape for each domain, and then offset reduced forces so that the new shape is the rest configuration. 4.4 Fracture Figure 4.9: Real-time Fracture: The tea bush (Camellia Sinensis) leaves are shaken from the bush by the user-applied force in real time. Left: Before fracture, with applied user force indi- cated. Middle: during fracture. Right: after fracture. Note that the leaves close to the user force location fractured in greater numbers than elsewhere on the model because they underwent higher accelerations. 6 msec of simulation per graphical frame. In nature, plants often fracture at the joints between its parts. For example, leaves or fruits detach from branches (see Figure 4.9), or branches crack away from the main stem. Such fracture with pre-specified patterns is useful in interactive applications because it is control- lable and artist-directable [74]. We support such fracture by monitoring the (reduced) interface forces between adjacent domains. If the L 2 -norm squared of reduced interface force vector ex- 4.5 Motion of Leaves 43 ceeds a user-adjustable threshold, we fracture the entire subtree from the main structure. The subtree (in many applications a single domain) then undergoes a ballistic trajectory under grav- ity, e.g., peach tree fruits land on the ground (Figure 1.2). Such fracture is computationally extremely inexpensive as the L 2 -norm test can be performed in nanoseconds. Because each part has its own rendering mesh independent of all the other parts, the objects are automati- cally free of holes after separation; there is no need for any re-meshing. The fracture events could also be scripted / keyframed, e.g., to simulate trees undergoing an explosion. Our do- mains can only fracture at the interfaces to other domains. For more detailed fracture, domains can be divided during the pre-process, e.g., the trunk can be pre-cut into two pieces. 4.5 Motion of Leaves Figure 4.10: Leaf Motion: A leaf is swinging in a randomized wind force field. Leaves are usually (but not always) considered as rigid domains in our method, especially in the tree simulations (Figure 1.1, Figure 1.2). It is possible to make leaves flexible by treating them as every other domain. (Figure 4.11, Figure 3.4) A rigid leaf has no reduced degrees of freedom so they can do only rigid transformation along with the branch it attaches to. This makes it possible to produce very plausible simulation results. However, to make the animation even more realistic, we can add motion to the leaves not by actually simulating them but by kinematically deforming them before the rendering stage. Note that this is an alternative approach to treating leaves as flexible domains.(Figure 4.10) To do so, we first compute a 3- dimensional linear reduced basis U ∈R 3n× 3 for a leaf. Here, n is the number of vertices of the leaf mesh. Since the leaves are instanced, the reduced basis is only computed for every representative leaf. Then we introduce a space-time randomized wind field whose strength is adjustable by the user. At each timestep, we query the wind field with the current time instant 4.6 Results 44 and the position of a sample point (the anchor point of the leaf), yielding a 3-vector f which is the wind force vector. We assemble the displacement of the leaf using: u( f)= U f (4.1) Sometimes the deformation of a leaf could be large, causing visible artifact since we are us- ing linear modes. We fix this problem by introducing the large deformation correction using modal derivatives [8]. Equation 4.1 essentially gives a kinematic term that corrects linearized deformations to large deformations. u( f)= U f+ 1 2 3 ∑ i=1 3 ∑ j=1 Φ i j f i f j (4.2) whereΦ i j is the modal derivative andΦ i j =Φ ji 4.6 Results We pre-processed over 100 plants; we provide a selected subset of 32 pre-processed and sim- ulated models in Figure 4.12. It can be seen that the performance is interactive even for very complex plants. Plants of small to moderate complexity are fast enough to be easily used in real-time systems such as computer games. Inverse kinematics and fracture were demonstrated in Figures 4.8 and 4.9, respectively. Plant motion resulting from user (mouse) forces, followed by free vibration, is demonstrated in Figure 4.13. Several of our plants are animated deforming in the wind. We refer the reader to Chap- ter 6 for more details about the wind. Our wind consists of two components: a wind with a (tunable) constant direction and magnitude, and a randomized wind, implemented as a 4D space-time Perlin noise [77], with standard parameters: number of frequencies and persistence (how quickly high frequencies decay). All the wind parameters are easily adjustable at run- time without any precomputation. Our method can model the entire spectrum of winds from gentle breezes, moderate winds (Figure 1.2) to hurricanes (Figure 4.14). Stochastic noise is 4.6 Results 45 Figure 4.12: Representative subset of simulated models: including flowers, bushes, broadleaf and conifer trees. Input meshes are from Xfrog, SpeedTree and 3dmolier (Turbosquid) model libraries. Models were deformed either by pulling on vertices or using a randomized wind. Each plant reports #triangles, #domains, #flexible domains, total # of DOFs, the simulation frame rate and memory. The simulation frame rate includes all computation to produce the next graphical frame, except the rendering itself. The pre-processing times range from a minute for simple models to 20 minutes for the most complex trees. Intel Xeon 2x8 cores 2.9 GHz CPU, 32GB RAM. GeForce GTX 680, 2GB RAM. Models from Table 4.1 are shown in (1-indexed) (rows, columns) (1,4), (2,1), (8,3), (5,3), respectively. 4.6 Results 46 Figure 4.13: Space-time instantaneous (localized) external force followed by free vibra- tion.Secondary motion in the smaller branches due to motion of main trunk is clearly visible. Simulation time: 5 msec per graphical frame. often used to animate trees by directly driving the deformations [73]. We employ Perlin noise to create spatially varying and controllably turbulent wind forces, causing large deformations and secondary motion of branches due to the motion of parent branches. We use Perlin noise directly, but other noise generators such as the 1/ f β noise [73] or even Navier-Stokes equa- tions [1] could be used instead. In order to apply the wind, we must sample it at properly selected locations on the model. One could evaluate the wind at every plant vertex, but doing so would require many wind evaluations and subspace force projections. Another alternative that we considered but did not pursue due to prohibitive cost was to sample the wind on a regular grid, and then seek the nearest vertex to each grid point. Instead, we select a repre- sentative set of volumetric mesh vertices on each branch domain, and then sample the wind at those locations (see Figure 4.14, top-left). Each wind sample is scaled with the volume of each branch. Our system also supports plants loaded by gravity (Figure 4.15). Because gravity acts in a constant direction, it needs to be rotated into the frame of reference of each domain, f ext i =U T i R i f, where U i and R i are the matrix specifying the low-dimensional space, and world- 4.6 Results 47 Figure 4.14: Palm tree (Cocos Nucifera) in strong randomized wind: Top-Left: the locations (in red) where the wind is sampled. Other images: selected animation frames. Simulation time: 4 msec per graphical frame. In this demo, palm leaves are not simulated, but are skinned to the palm branches. Each leaf is skinned entirely to one branch, with each vertex copying the displacement (in local branch frame) of the closest branch vertex in the undeformed configuration. Such skinning can greatly enrich the plant visual appearance at a minimal computational cost. Figure 4.15: Gravity space:Left: undeformed. Right: under gravity. 4.6 Results 48 coordinate rotation of domain i, respectively, and f =[0,− g, 0, 0,− g, 0,...0,− g, 0] T is the gravity vector. Because f is constant, f ext i can be efficiently precomputed using the sandwich transform [48], by evaluating f ext i for R i = e k e T ℓ , where e k is the k-th standard basis vector in R 3 , for all k,ℓ= 1,2,3. At runtime, for each non-rigid domain, one then merely has to multiply a constant pre-computed r i × 9 matrix (r i = #reduced DOFs of domain i) with the 9-vector of the entries of R i , imposing a negligible overhead. It would be possible to pre-load plants so that the input configuration is also the rest configuration under gravity [97]. Figure 4.11: Visualization of domains in the oak tree: In this example, all leaves are flexible domains, each of which has 8 reduced DOFs. Comparison to full simulation: In practice, vegeta- tion is often simulated by building a simulation mesh for the entire model, then employing a deformable object simulator to timestep the plant dynamics for- ward in time [4, 59]. In Figure 4.16, we compare our method to a geometrically nonlinear full sim- ulation [16]. For this experiment, we generated a global tetrahedral mesh for the shefflera plant (Shef- flera Actinophylla), and then manually subdivided it into domains. This process took 8 hours of work, whereas our pre-processing pipeline takes less than 10 minutes. In order to make the motion more natural, we made the stem 2x stiffer than the rest of the mesh. We then simulated the shefflera under an identical force load, material and simulation properties. In this exam- ple, our method (including time for u= Uq displace- ment computation) is 23x faster than the full simula- tion. Visually, the two motions differ slightly in frequency, but appear qualitatively similar. If a close frequency much is desired, it is possible to automatically scale the frequency spectrum of each domain so that the lowest frequency matches some externally prescribed frequency, such as frequency from full simulation, or real measurements of plants [43]. We illustrate this 4.6 Results 49 concept with a 1D harmonic oscillator, m ¨ x+ d ˙ x+ kx= f, whose natural angular frequency (without damping) isω = p k/m. Suppose a different frequencyω ′ ̸=ω is desired. The mod- ified mass and stiffness must therefore satisfy α 2 km ′ = k ′ m, where α =ω ′ /ω. If we impose an additional condition that the new oscillator attains the same maximum amplitude after a fixed initial impulse, we obtain k ′ = αk and m ′ = m/α. With plants, we perform the same scaling, to match the lowest vibration frequency, separately for each domain. We obtain the desired frequency for each branch by creating a volumetric mesh for the entire subtree rooted at that branch, and computing (unreduced) mass and stiffness matrices M and K. The desired angular velocity is then the square root of the lowest eigenvalue of Mx=λKx, which we find using a sparse eigensolver [54]. We scale the reduced mass matrix by 1/α and internal elas- tic forces by α. The result is shown in Figure 4.16. We note that, instead of matching the lowest frequencies to full FEM simulation, one could instead match them to real-world plant observations. Figure 4.16: Comparison to full simulation. 91 flexible domains, 1183 modes, 11,396 trian- gles. All three simulations visually look similar. The unscaled reduced simulation has a higher natural frequency, for two reasons: (1) reduced systems lack degrees of freedom and are typically somewhat stiffer than unreduced systems with same parameters, and (2) interface lumping [Barbiˇ c and Zhao 2011] increases frequency, much like a pendulum with a shorter length oscillates faster. After scaling, the frequency and amplitude match full simulation closely. 4.6 Results 50 Rendering: We render the plants interactively using OpenGL. Billboards are in practice of- ten partially transparent, with alpha values continuously ranging from zero to one. Therefore, one-pass alpha-testing results in noticeable aliasing, for example, at the leaf edges. We use two-pass rendering where we first render all the geometry with alpha-testing enabled, i.e., fragments where the alpha value is strictly less than 1.0 are discarded. We then make the depth buffer read-only, disable alpha-testing, enable alpha blending, and re-render all domains that use transparent texture maps. Note that in many models transparent textures represent a large fraction of the geometry, e.g., leaves, conifer needles, twigs. Our simulation is substantially faster than rendering; more optimized rendering pipelines could be designed [90]. Before our models can be rendered, the displacements (in local domain frame of reference) of all the mesh vertices must be computed via the modal equation u i = U i q i [40] and then their world coordinate positions constructed using the domain’s current position vector and rotation ma- trix. All our reported simulation timings include the time necessary to compute u i and world coordinate triangle mesh vertex positions. We compute u i on the CPU, but GPU implemen- tations [40] would be readily possible. We note that matrix U i here contains the modes that were already interpolated (during pre-process) from the volumetric mesh to the plant triangu- lar geometry, using barycentric interpolation. We found that such a strategy is usually faster than interpolating volumetric mesh displacements to the triangle mesh at runtime; the tradeoff depends on how finely the triangle meshes are tessellated. Offline renderings were performed using the Yafaray ray tracer. Chapter 5 Procedural Materials and Asynchronous Parallel Timestepping Given a large plant, how to determine material properties of each branch, twig or leaf? In real trees, branches undergo lignification, which greatly affects mechanical behavior of trees. In the first part of the chapter, we present a novel procedural method to automatically set sim- ulation parameters for complex plants. These parameters include material stiffness, damping coefficients, and timestep. Such tuning is quite challenging due to the complexity of realistic tree models. In practice, it is often too complicated to do it manually in order to generate vivid and plausible plant motion. The simulation results, achieved by applying automatically tuned parameters to our reduced multi-domain simulator, demonstrate the quality and efficiency of our tuning algorithm. Note that this algorithm could be applied not only to our architecture, but also to any other FEM simulation. In the second part of this chapter, we introduce an asynchronous, parallel timestepping algo- rithm to improve our reduced multi-domain dynamics simulator presented in previous chapter. Since the natural frequency varies in domains, we integrate each domain at its own pace, as opposed to using a fixed timestep for the entire object. This avoids the limitation of a syn- chronous integrator, in which a small, fixed time step needs to be carefully chosen in order to avoid numerical instabilities due to the high stiffness variations encountered at the domains. The total running time to produce animations of a given length (same total physical time) was 5.1 Procedural Stiffness 52 therefore shorter than the synchronous method. By exploiting the fact that no loop exists in the tree topology, the simulation is further accelerated by timestepping all domains in parallel. 5.1 Procedural Stiffness We now present our approach to set the material property of each domain so that the result motion looks real, and natural. In nature, objects vibrate at different frequencies when their size is changed. Longer and thicker objects tend to vibrate slower. If we keep the material properties constant and uniform everywhere on the plant, smaller branches would then have much higher natural frequencies than the bigger ones. This will cause, for example, such small branches to be severely overdamped and behave mostly rigidly. However, this does not match the phenomenon observed in the real world where smaller branches usually tend to behave more lively. For example, twigs are swinging heavily in the wind while the trunk may not shake too much. The reason is that branches have different material stiffness. Small branches are usually younger and therefore feels softer than the big ones which have been largely lignified. Therefore small branches in nature have lower vibration rates. To create realistic results, one has to set the stiffness of all branches correctly. A possible way to do this is to adjust the stiffness of each branch manually. The more dense the object is, the slower it will vibrate, and the lower its frequency will be. Such a technique, however, requires a lot of manual work and is by no means intuitive to the user. In this section we design a model called “aging mode” which automatically sets the stiffness of each domain based on the geometry of a domain and the topological structure of the entire plant. We observe that the age of a branch often relates to its length. In most cases, the older the branch is, the longer it will be. Our idea is to first estimate the age based on the geometry and then set the stiffness according to the age. To do so, we design the following scalar function to model the stiffness as a function of age. E(a)= c 1 · a D (5.1) 5.1 Procedural Stiffness 53 where a is the age of the branch and c 1 (c 1 ≥ 0) and D (D≥ 0) are two real coefficients. The exponent D is input by the user from UI. It controls how aggressive the aging mode is. Higher values will make small branches younger and therefore more lively. Stiffness is proportional to the square of the lowest natural vibration frequency. Therefore we model the final branch frequencyν as ν =ν geom · p E(a) (5.2) whereν geom is the initial frequency under default, uniform materials, precomputed in the pre- processor using linear model analysis. By inserting 5.1 into 5.2, we obtain ν = c· ν geom · a D/2 , (5.3) where c= √ c 1 . In the following, we explain how to model the age of the branch. Letℓ denote the length of the branch. We postulate that the age is linearly proportional to lengthℓ. Since ν geom is inversely proportional toℓ, we have: a= const· ℓ= c 2 ν geom , (5.4) where c 2 ≥ 0 is a constant. Sometimes a branch is shorter than its descendents. That means it has a higher initial frequency than some branches in its sub-tree. This often happens, for example, in “shrub models”. In this case, Equation 5.4 implies that this branch is “younger” than some of its descendents which is a contradiction. To fix this problem, we enforce the condition that a branch should be always order than all of its descendent. We extend 5.4 to a= max{ c 2 ν geom , a oldest_descendent } (5.5) where a oldest_descendent is the age of the oldest descendent of that branch. Our material adjusting procedure consists of two steps. In the first step, we traverse hierarchy from leaf nodes to the root (bottom-up) and compute the age property for every domain based on Equation 5.5. In the second step, we set the stiffness by assigning vibration frequency of 5.2 Procedural Timestep 54 each branch based on Equation 5.3. We set the coefficient c constant within each hierarchy level. User sets c for each level on the UI. This is based on the assumption that domains with the same depth belong to the same generation and have roughly the same stiffness property. We assign decreasing c values to increasingly deep levels from root to leaves. A good way to tune a tree is to set the first few levels stiff, i.e., 3.0, 5.0 or even 10.0. Then, make deeper levels softer, i.e., 1.0. For the aging exponent D, a value of 2.0 is a reasonable choice; it causes all branches to have equal frequency. A value of 0.0 is equivalent to a constant stiffness for all the branches, i.e. purely geometric models. From the simulation results, for D= 2 (or similar) we clearly observed much richer and more lively secondary motion in smaller branches due to the motion of parent branches. And the main body of the plant, like the trunk, looks very stiff preventing it from rubber-like deformations even under large external forces. The overall motion looks much more natural and real then when using a constant and uniform material property everywhere in the tree. 5.2 Procedural Timestep Here our goal is to compute a separate timestep for each domain so the simulation could be advanced in each local frame without any numerical instabilities. To achieve this, we first solve a mass-scaled generalized eigen problem [8] to compute the natural frequencies of the domain. The generalized eigenproblem has the form, Kψ =λMψ where the K is the reduced tangent stiffness matrix, evaluated at the origin. Here, M is the reduced mass matrix (Equation 3.17). Note that both K and M are r× r symmetric, positive-definite matrices, where r is the number of degrees of freedom of the reduced space. If the r eigenvalues are sorted in ascending order, the smallest eigenvalueλ 0 > 0 corresponds to the lowest natural frequency, also known as the main frequency, of the object. We pick thisλ 0 to compute the corresponding oscillatory period as λ =ω 2 where ω denotes the natural frequency in radians. The oscillatory period is then given by T = 2π p λ 0 . (5.6) 5.3 Procedural Damping 55 At this point, our system requires the user to provide a positive integer, called quality control coefficient , denoted as Q. The simulation timestep is set as ∆t 2 N . where n satisfies that ∆t 2 N < T Q < ∆t 2 N− 1 , (5.7) and ∆t is the graphical time step. We usually set ∆t = 1/30 sec or ∆t = 1/24 sec, which conforms to the standard frame rate in modern TV and digital cinema. The user is able to adjust this value on the UI. We use power of 2 (2 N ) so that we can perform parallel timestepping efficiently (next section). The quality control coefficient Q is the only parameter that needs to be set by the user. By default, we set Q= 100 which is, in our experience, a very good value for most simulations. User can easily adjust it to favor either the simulation speed (lower Q) or accuracy (higher Q). Our simulator is quite stable even for small values of Q where a lot of artificial damping is introduced into the system which helps to quickly dissipate the energy. We only observed explosions in some very rare and extreme cases. By increasing Q, these numerical instabilities were soon fixed and disappeared. 5.3 Procedural Damping We adopt proportional (Rayleigh) damping in our system which has the form C=αM+βK where M is the reduced mass matrix, K is the reduced stiffness matrix, α and β are damping coefficients. Since domains have different stiffness and are simulated separately, each of them should be assigned different damping coefficients α andβ in order to achieve better simulation results. To do so, we have designed a procedural method to tune these damping coefficients automatically. We first employ the modal damping factor [40]: ξ = 1 2 α ω +βω (5.8) 5.3 Procedural Damping 56 whereω is the undamped, lowest natural frequency of vibration (in radians) computed using: ω = 2πν, (5.9) where ν is the lowest natural frequency of the domain. We model ξ as a function of the frequency of each domain. This is done by modelingξ as a constant function ofν forν≤ ν low andν≥ ν high and a straight line in between. Formally, this piecewise function is as follows: ξ(ν)=        ν low :ν≤ ν low ξ high − ξ low ν high − ν low (ν− ν low )+ξ low :ν low <ν<ν high ν high :ν≥ ν high (5.10) whereξ high ,ξ low ,ν high , andν low are set by the user on the UI. Given the frequency of a domain, ξ is then computed using equation 5.10. However, knowingξ does not uniquely define α and β since we have two unknowns but only one equation 5.8 so far. At this point we introduce another dimensionless factor γ which models the liveliness of the high frequency components. When the natural vibration frequency becomes twice as large as the lowest one, we let the modal damping factorξ be scaled by 1/γ, that is, 1 2 α 2ω + 2βω = ξ γ . (5.11) The damping coefficients thus can be computed using Equations 5.11 and 5.8: α = 4ξω 3 (2− 1 r ), (5.12) β = 2ξ 3ω ( 2 r − 1). (5.13) Generally we wish the damping forces to penalize higher frequency motions more, which meansξ should be a monotonically increasing function for frequencies of each domain. It can be derived from Equation 5.8 that ξ reaches a global minimum when ω ∗ = q α β . We impose a condition that whenξ reaches the minimum, the frequencyω ∗ is no greater than the lowest 5.3 Procedural Damping 57 frequency of the domain (Figure 5.1), that is, ω ∗ = r α β ≤ ω (5.14) By inserting Equation 5.12 and 5.13 into 5.14, we obtain the condition γ ≤ 0.8. Since α and β are non-negative values, we get γ ∈[0.5,2]. So eventually we have γ ∈[0.5,0.8]. On the UI, the user inputs a value γ user within the range [0.0,1.0] which is mapped to γ by γ = 0.5+ 0.3γ user . By default,γ user is set to 0.0 which leads toα = 0.0. This means that only stiffness damping β is injected into the system which is typically sufficient to stabilize the simulation. Increasingγ user will add more high frequency motions to the simulation. Figure 5.1: Avoidξ ’s minimum The minimum point ofξ is marked in yellow. The lowest natural frequency (in radians) of the object is represented by a red dash line. We impose thatω ∗ ≤ ω. 5.4 Asynchronous Parallel Timestepping 58 5.4 Asynchronous Parallel Timestepping Synchronous time integration evolves the entire plant forward in a fixed step from one time instant to the next. It is simple and straightforward. However, the achievable simulation speed is dictated by the stiffest domain which often requires a very small timestep. Such timestep, however, is unnecessary for the domains whose vibration frequency is much smaller. These domains could be time-evolved more efficiently using a much bigger timestep. Therefore, we design an asynchronous integrator which advances different domains at different rates. 5.4.1 Subdividing the Timestep As presented in Section 5.2, the user is required to set the graphics timestep∆t, and the quality control coefficient Q . The physics timestep will be divided into 2 N sub-timesteps per domain according to Equation 5.7. This sub-timestep is the actual simulation step. Our first attempt was to modify Algorithm 1 so that we would do semi-implicit Newmark integrations N times. It was performed by incorporating Line 9 in Algorithm 1 into a for loop from 0 to 2 N− 1 . Note that we did not even need to use powers of 2 at this point. We only needed a factor N to satisfy ∆t N ≤ T Q . This approach, however, did not work very well in practice. Because the local frame (polar decomposition), relative and absolute kinematics were updated only once before the in- tegrations. Vibrations and instabilities could occur for large deformations and small values of Q, as the integrator was trying to move forward with too many steps using obsolete kinematic and dynamic information. To fix this problem, we set the subdivision factor to be a power of 2 which satisfies Equation 5.7. This guarantees that we can update domain’s kinematics and dynamics uniformly in time. Note that if domain is rigid, we set N = 0 if it is also the root, otherwise, we set N equals to its parent’s N. If a sub-tree is found rigid, namely, all domains within are rigid, we set N = 0 for all domains in the subtree. This always happens when a large number of leaves recognize a twig as their parent. In such cases, this saves a significant amount of computation time. 5.4 Asynchronous Parallel Timestepping 59 We update the absolute and relative kinematics once everyℓ timesteps, whereℓ is also a power of 2 and varies with each domain. We will explain how to compute ℓ later in this section. Assume the domain being simulated is domain i. When the integrator is at the beginning of eachℓ timestep sequence, it will first update (i) absolute kinematics, (ii) system forces, and interface forces, then (iii) doℓ semi-implicit Newmark integrations forward, (iv) compute local frame and relative kinematics for every child domain, and finally (v) save updated absolute kinematics v i ,a i ,ω i ,α i to domain i’s buffer, and relative kinematics x i j ,A i j ,R i j ,v i j ,a i j ,ω i j ,α i j (see Section 3.2) to the buffer of its child domain j. The buffer size for absolute kinematics of domain i is 2 N i /ℓ i while the buffer size for relative kinematics of domain j is 2 N j /ℓ j . Using the buffer prevents the parent from waiting for the child in the parallel algorithm presented later in this chapter since the child needs kinematic information from the parent. If domain is rigid, steps (ii),(iii) are skipped. Now we describe how to determineℓ. We assume a general case where domain j has a parent domain i and n child domains, k 0 ,k 1 ,...,k n− 1 . Thenℓ j is computed as: ℓ j = min{ 2 N j − N i , j 2 N j − N k 0 k , j 2 N j − N k 1 k ,..., j 2 N j − N k n− 1 k } (5.15) where N i ,N j ,N k 0 ,N k 1 ,...,N k n− 1 are the exponents for these domains. If domain j is the root (no parent) or a leaf (no children), we simply remove the related terms in Equation 5.15. Figure 5.2 illustrates a concrete example of our timestepping algorithm working on a tree. 5.4.2 Parallel Simulation Since no loop exists in the tree structure, it is possible to simulate domains in parallel. The fundamental idea is as follows. We setup simulation tasks and push them in a work queue. The task starts when the scheduler pops it from the queue. It happens when a CPU thread (core) becomes available. Note that the work queue is a global, critical resource and the “push/pop” operations must be serialized. Once created, all CPU threads are kept alive, executing or waiting the tasks until all the work is done. Each domain has 2 N simulation steps. Based on the content of the tasks, we have designed and implemented two scheduling schemes: 5.4 Asynchronous Parallel Timestepping 60 Figure 5.2: Asynchronous Timestepping: Each independent rectangle represents a simulation step of one domain. Rectangles have the same width indicating the graphics timestep (e.g., 1/30 sec) is fixed for every domain. Every rectangle is subdivided into 2 N grids where 2 N is the subdivi- sion factor of the domain. Each grid represents a physics timestep for a domain. An arrow pointing from domain i at time instant t 0 to domain j at time instant t 1 means that the domain j requires domain i’s kinematic information at time instant t 0 to update its own kinematics at time instant t 1 . The update frequencyℓ is shown beside each rectangle. Scheme 1: Atomic task (fundamental amount of work by each thread) is to simulate a do- main 2 N physical steps (N varies with each domain, as discussed in previous section). This means simulating a single domain occurs exclusively on one thread. Data dependency hap- pens when a domain needs its parent’s kinematic information. This will not be an issue since we adopt breadth first search for the tree traversal. Domain is fully simulated prior to any of its children. In this process, domain assumes the parent has provided the necessary data to the proper internal storage location, and then reads from it. When the task is completed, all the child domains, which can be timestepped in parallel, will be pushed into the work queue. Algorithm stops after all domains are simulated. Scheme 2: Atomic task is to simulate a domainℓ physical steps (ℓ≤ 2 N , see Section 5.4.1). The collection of tasks is represented by a directed acyclic graph (DAG) (Figure 5.3) with a node for each task and an edge for each constraint (data dependency). The DAG gives an order when certain tasks must be performed earlier than others. Compared to the first scheme, there is one more data dependency here, that is, the task to integrate domain from step k (k≥ 1) cannot start if step k− 1 has not been computed. In order to generate a valid task sequence, we 5.4 Asynchronous Parallel Timestepping 61 Figure 5.3: Directed acyclic graph: Each node is a simulation task while each directed edge is a constraint representing the data dependency. employ the standard topological sort algorithm. Traversal starts from simulating root domain at timestep 0. Once the job is done, program writes all intermediate results that may be needed by the following tasks internally to the buffer. A new task is triggered immediately when all its constraints (there are at most two) are satisfied. Comparison The second scheme provides a more compact schedule plan than the first one. Theoretically, it means the simulation is expected to be faster. However in practice, we found it is slower than the first scheme in most cases. The reason is buffers storing the kinematic information may cause read-write conflict between child and parent domains. To prevent this conflict, additional locking is needed which imposes additional complexity and overhead to the algorithm. When the number of tasks is large, the lock/unlock operation happens fre- quently, causing significant decrease in performance. As opposed to the second scheme, the first scheme does not require locking. By default, our simulator uses the first scheme. User is allowed to switch to the second scheme dynamically during the run time. The system perfor- mance has been compared and illustrated in Table 5.1. 5.4 Asynchronous Parallel Timestepping 62 Example triangles domains flexible DOFs FPS (SIGGRAPH2013) Quality FPS (Scheme 1) FPS (Scheme 2) Tea bush 12,730 301 37 528 166Hz 1 1,604Hz 1,428Hz Conifer 5,587 644 43 360 316Hz 50 640Hz 586Hz BroadLeaf HD 38,272 7,003 420 3,613 41Hz 100 52Hz 49Hz Fir 81,228 9,077 588 4,858 25Hz 150 50Hz 40Hz White oak tree 2,360,868 120,871 871 9,520 1Hz 50 6Hz 5Hz Table 5.1: Simulation statistics for #triangles, #total domains, # flexible domains, total # of reduced DOFs, frame rate (fps) from our SIGGRAPH 2013 paper [106], quality control coefficient, frame rate (fps) for Scheme 1, domain-based multicore algorithm, frame rate (fps), for Scheme 2, DAG-based multicore algorithm Machine specs: Intel Xeon 2× 8 cores 2.9 GHz, 32 GB RAM, Nvidia GeForce GTX 680, 2GB RAM. Chapter 6 Randomized Wind Force Field Wind plays a very important role for adding natural-looking variety to plant animations. This chapter describes how we create a randomized wind force field used in our simulator. Our im- plementation is based on the improved Perlin noise algorithm presented by [77]. We also refer the reader to [89] as good supplementary materials. We begin by describing the requirements for a correct and robust randomized wind field. Next, we provide implementation details on the 4D Perlin noise algorithm. Then we extend technique to construct a randomized wind. Finally, we demonstrate how to apply the wind to the plants in our system. 6.1 Features of a Randomized Wind Force Field Wind is an important building block in plant simulation. A good wind force field incorporates (but is not limited to) the following features [89]: • It is not spatially or temporally uniform, namely, it has a controlled way of adding ran- domness to the output. • It produces a repeatable pseudorandom force vector for every input pair, position and time instant. • It is smooth which means it has band-limited spatial and temporal frequency. 6.2 Improved Perlin Noise Algorithm in 4D 64 • It does not show obvious repeating patterns. • Its spatial frequency is invariant under translation. The improved Perlin noise [77] satisfies these requirements, and we adopt it in our work. 6.2 Improved Perlin Noise Algorithm in 4D Figure 6.1: Visualizing the Perlin wind field: The bounding box of the wind field is shown in red. The spacebox of the wind field has been subdivided at a res- olution of 4× 4× 4. We sample and vi- sualize the Perlin wind (blue line seg- ment) at every grid point (shown in green, 5× 5× 5). To improve visualization, the noise contains a single frequency compo- nent. The standard improved Perlin noise algorithm takes a 3D position as input. We extend it to incorporate time as the fourth dimension. In general, user inputs a real quadruple (x,y,z,t) where (x,y,z) is the 3D position and t is the time instant. The noise at this 4D position is obtained by first computing a pseudorandom gradi- ent at each of the sixteen nearest quadruple nodes (the “corners”) on the integer 4D lattice, and then perform- ing spline interpolation (see [77]). Let (x ′ ,y ′ ,z ′ ,t ′ ) denote the coordinate of one corner where x ′ =⌊x⌋ or⌊x⌋+ 1, y ′ =⌊y⌋ or⌊y⌋+ 1, z ′ =⌊z⌋ or⌊z⌋+ 1 and t ′ =⌊t⌋ or⌊t⌋+ 1. The first step is to compute a pseudorandom gradient for(x ′ ,y ′ ,z ′ ,t ′ ). We create and store a hash table which contains a random permutation of integers from 0 to 255. This hash table is indexed by the x ′ coordinate first. The value at po- sition x ′ is added to y ′ coordinate and the sum is used to index the hash table again. This process is repeated until t ′ is incorporated. The result is a pseudorandom integer for the corner(x ′ ,y ′ ,z ′ ,t ′ ). This pseudorandom integer is used to look up into a table of 4D gradient vectors. The content of this 4D gradient table is precomputed in the initial stage. We refer the reader to [89] for more details. In the second step, we obtain a noise value for the corner(x ′ ,y ′ ,z ′ ,t ′ ) by computing a dot product between the gradient vector and the fractional 6.3 Computing the Randomized Wind 65 position (x− x ′ , y− y ′ , z− z ′ , t− t ′ ). Finally we achieve the output value by quadrilinearly interpolating the sixteen noise values using blending factors w(x− ⌊x⌋), w(y− ⌊y⌋), w(z− ⌊z⌋) and w(x− ⌊z⌋), where w(t)= 6t 5 − 15t 4 + 10t 3 . 6.3 Computing the Randomized Wind Our wind consists of a constant wind (a 3-vector) with the given direction and magnitude, and a randomized spacetime Perlin noise (a 3-vector). Such a randomized noise is able to model strong, turbulent winds that cause our plants to deform to large deformations. The noise is made up of several components at different frequencies. If they are sorted in ascending order by frequency, the contribution of each component follows an exponential decay function, that is, g(i)=(1− d) i where d is the decay rate within the range[0,1], and i is the frequency index and g(i) is the contribution of the i th frequency. Therefore the input to our algorithm includes the direction and the magnitude of the constant wind, the position of the sample points(x,y,z), the time instant t, the initial amplitude of the randomized wind A, the number of frequencies, and the decay rate d. We now explain how to compute each component of the randomized wind. This is done by concatenating the sample spatial position and time instant, scale it with the current frequency, and provide it to the Perlin noise generator (Section 6.2) to generate an initial noise. We multiply the initial noise with the contribution g(i) and the initial amplitude A to obtain the final noise. We add this final noise to the constant wind to form the randomized wind. We repeat this process for the other 2 components of the wind. 6.4 Applying Wind Forces to the Plant We compute the axis-aligned bounding box of the scene and scale it by a ratio (by default, ratio = 1.2). The new bounding box defines the space of the randomized wind field. To visualize such a vector field, we first divide the space into cubic voxels based on user-defined resolution. Then we sample the wind force at every grid vertex and render it on the screen. 6.4 Applying Wind Forces to the Plant 66 Figure 6.2: Applying Perlin wind to the plant: The tea bush (Camellia Sinensis) is placed in a Perlin wind field. (a) The tea bush is swinging in the wind. (b)(c) Visualization of the wind forces on sample points of the tea bush. The sample points are shown in red. The wind forces applied on sample points are shown in green. For better visualization quality, only one point is sampled per domain (number of samples: s= 1) (d) Visualization of the Perlin wind field. The noise contains 3 frequency components, and the decay rate is 0.05. 6.4 Applying Wind Forces to the Plant 67 (Figure 6.1) In order to apply the wind to the simulation, we sample vertices on the volumetric mesh where the wind forces will be applied. A naive idea is to use all the vertices on the surface of the volumetric mesh. This is not practical since Perlin wind evaluations and reduced force projections are computationally expensive and the number of sample points is too large. The solution adopted in our system is as follows. The user supplies the desired number of samples s≤ 1 per domain. We use a constant number of samples for all the branches (typically s= 5), but the number could be scaled with domain size, e.g., s= 1 for leaves. The sampled vertices for a domain are then determined automatically, as follows. We first build a mesh graph for our voxel domain mesh, where nodes are volumetric mesh vertices and nodes are connected if they are adjacent mesh vertices. We then use Dijkstra’s algorithm to compute the minimum graph distance of each volumetric mesh vertex to the set of fixed vertices for this domain. Let D denote the maximum distance. Then, for each i= 1,...,s, we select any vertex with distance⌊iD/(s+ 1)⌋ as our sample. Because plant branches are long and slender, and the running time of Dijkstra’s algorithm is linear in the number of vertices of each domain, such a strategy produces well-distributed samples in negligible time. In order to model the fact that large branches are more exposed to wind, the wind magnitude for a sample point is set to mass/s where mass is the mass of the volumetric mesh of each domain. To apply wind forces to the plant, we compute the randomized wind force vector on every sample point on every domain and accumulate it in our simulator (Figure 6.2). In Figure 6.3, we demonstrate the motion of a grand fir subject to the Perlin wind. 6.4 Applying Wind Forces to the Plant 68 Figure 6.3: Grand fir (Abies Grandis) in Perlin wind: Frame numbers are provided. Chapter 7 System Overview In this chapter, we provide an overview of our system which consists of two applications. The first one is the botanical preprocessor and the second one is the botanical simulator. We introduce the graphical user interface (GUI) for each program separately and describe how to tune the parameters on the GUI. 7.1 Botanical Preprocessor 7.1.1 Command Line User can use the command-line (no GUI) if the hierarchy is already known. Pass the hierarchy via -h. User has to also pass -a to avoid the GUI. In this way, user can pre-process the tree completely in the shell, without any user interaction. The hierarchy file has a simple text format where each line lists a domain and its parent. Root is assigned -1 as its parent. You can save a hierarchy from the GUI by pressing “Save”. 7.1.2 Preprocessor User needs to use the GUI if the hierarchy is not known or if it is not connected or contains cycles. When the program is launched, a panel is shown (Figure 7.1). 7.1 Botanical Preprocessor 70 Figure 7.1: Plant preprocessor Left: Botanical object display window. Right: Preprocessor UI panel Selecting fronds and leaves This makes it possible to separate the domains into branches, fronds and leaves. If all domains are branches, user can simply press the “Frond & Leaf Selection Done” button. Otherwise, select some domains, and press ‘f’ to tag them as fronds or ‘l’ as leaves (Figure 7.2). The user can delete domains by pressing the ‘del’ key. Selecting root Select the root domain and press “Set Current Group As Root (F5)” button or ‘F5’ (Figure 7.3). Building hierarchy This button builds the hierarchy, by performing collision detection on the triangle meshes. If more than one connected component is detected, and/or there are loops, additional steps are necessary (below). Connecting the domains The GUI will display the connected components. Scroll via ‘,’ and ‘.’. User can delete domains by pressing the ‘del’ key. To connect two disconnected components, select the first domain (LMB) from the first component, press ‘ c’. Then, select 7.1 Botanical Preprocessor 71 Figure 7.2: Select fronds and leaves (a) leaves are selected and highlighted in cyan, (b) After pressing ‘l’, leaves are tagged, (c) Fronds are selected and highlighted in cyan, (d) After pressing ‘f’, fronds are tagged. Figure 7.4: Connecting the domains: Select a domain from the first component and press ‘c’, the domain will be highlighted in blue. Then select another domain from the second component and press ‘c’ again. These two domains will be connected together. 7.1 Botanical Preprocessor 72 the second domain from the second component and press ‘c’ domain. The two domains are now connected in the domain graph. There is no undo for this. However, after pressing the first ‘ c’, user can press ‘C’ to cancel. Continue connecting until there is only one connected component. At any moment, user can press ‘Build Hierarchy’ and the hierarchy will be re-built. This way, user can see how many components were already fixed. User can press ‘ s’ to save the connectivity graph to the disk. It is saved to: output.reducedMultiDomainDeformationFactory.connectivityGraph<.suffix>. The “<.suffix>” is one of the following: 1. none: it is automatically saved by the system after doing collision detection to build the connectivity graph; 2. ‘.c + number’: c marks the program is in the resolving disconnected component stage, and the number indicates how many disconnected components are existing. For exam- ple: .c7; 3. ‘.b + number’: b marks the program is in the resolving loop stage, and the number indicates how many loops are remaining. For example: .b26; 4. ‘.allResolved’: it means all the ambiguities were resolved. This functionality allows the user to save the current work and continue preprocessing at a later time. To load the connectivity graph, use -G option and specify the filename when launching the preprocessor GUI via the command line. In this mode, a new button “Load” appears beside the “Build” button. User should go through all the previous steps described above and then click the “Build” to initialize the hierarchy for the first time. Then press the “ Load” button to load the connectivity graph and continue preprocessing from there. Resolving loops: If there are loops, they must be resolved manually using the provided in- terface. Our software uses an intelligent algorithm which minimizes the amount of user work required (see Section 4.1). As opposed to visiting every pair of colliding branches and asking the user if this is a genuine collision or not, the algorithm automatically selects a much smaller 7.1 Botanical Preprocessor 73 Figure 7.5: Resolving Loops: (a) Loop visualization. Branches, fronds and leaves are rendered. They can be hidden by pressing ‘1’, ‘2’, ‘3’. (b) Visualization of the loop while only branches are rendered. (c)-(e): cycle through the edges by pressing ‘{’ or ‘}’ and break the loop by pressing ‘d’. set of loops (call it S) in such a way that resolving those loops provably removes all loops. For each loop in S, the user needs to specify which edge in this loop is spurious. For example, if a branch accidentally collides with another branch, that would be a spurious edge. Note that an “edge” here refers to a pair of colliding branches. Nodes means branches while edges means colliding branches. (Figure 7.5) The spurious edge is selected using keys ‘{’ and ‘}’. You will see the branches in the loop in yellow, and the two branches whose joining edge is currently selected are shown in red and pink. If this is the edge to be removed, press ’d’. Otherwise use ‘{’ and ‘}’ to select the correct spurious edge, and then press ‘d’. User can also go to previous or next loop using ‘[’ and ‘]’, but this is typically not needed. Once all the loops are fixed, the program will allow the user to continue to the next stage (model reduction). User can undo with ‘z’. And the branches can be hidden/shown using ‘1’ which is useful to see just the loop. We summarize the keys for loop resolution: • ‘{’ and ‘}’: go to previous or next unresolved loop • ‘[’ and ‘]’: cycle through neighboring branch pairs within the current loop • ‘d’: break the cycle at the current branch pair (edge) (shown in red and pink) • ‘z’: undo the previous cycle resolving operation • ‘1’: hide or show the entire tree (keeping only the loop visible) • ‘s’: save the connectivity graph 7.1 Botanical Preprocessor 74 This functionality is also available on the UI panel. Note that once the hierarchy is a tree, user can save it to a file via the “ Save” button. It is saved to: output.reducedMultiDomainDeformationFactory.hierarchy. This hierarchy can then be loaded to the preprocessor via the ‘-h’ command-line option. (In this way, the user only has to perform connections and loop resolution once.) Note that the user can load the hierarchy and avoid UI as long as no domains needed to be deleted before the loading. If deletion is required, user has to use the UI. However, connectivity graph (if saved before) could still be loaded to save the preprocessing time. Also, once the hierarchy is built, the program will output the employed branches, fronds and leaves to files: out.<d>.obj whered is a digit and can be ‘0’, ‘1’ and ‘2’ corresponding to branches, fronds, and leaves. In this way, even when some domains were deleted, command line mode can still be used with -h. Figure 7.3: Select root domain: User first selects a domain (highlighted in cyan). And then presses ‘F5’ or the but- ton on the UI to set it as root. Model reduction After the domain graph has been connected and all loops resolved, user needs to press “Perform Model Reduction”. This will compute simulation meshes (via voxelization), determine the interfaces and perform model reduction as in [8]. At the end of this step, all the necessary simulation files will have been created. Now, you are ready to launch the runtime driver with the simulation config file as the the only input. Aging mode This option is accessible via the “Configure" rollout menu in the preprocessor. When the aging mode is on, branch stiffness is assigned pro- cedurally. It is enabled by default. If user turns it off, a tree has constant material stiffness everywhere. Note that the user can change this at runtime very easily (no additional preprocessing required). In the pre-processor, the user is only selecting whether the aging mode is enabled by default when 7.1 Botanical Preprocessor 75 launching the runtime driver. 7.1.3 Useful Tools In this section, we present some independent, useful tools along with the preprocessor. Prune a plant: The user is able to create plants obtained by keeping the first N levels of a plant, using the “prunePlant” utility. Given the input obj mesh and the hierarchy file, this will create the obj mesh that keeps the first N hierarchical levels. Create a forest: This utility helps the user to create a forest (see Section 4.2). It consists of following steps: 1. Pre-process the individual tree/plant (using the preprocessor) that will be replicated into a forest. 2. Create an instance description file (i.e., forest file, extension .forest). This text file looks like this: <num plants> x,y,z of plant 1 x,y,z of plant 2 ... 3. Use the utility “instancePlants” to create domain and parent list filenames for the for- est. Input parameters are the domain and parent list filenames created by the preproces- sor. Usually they are called “base-<obj name>-domainList.txt” and “base-<obj name>parentList.txt". For example: ./instancePlants base-myBigTree-domainList.txt\ base-myBigTree-parentList.txt myBigForest.forest\ base-myBigTree-domainList.forest.txt\ base-myBigTree-parentList.forest.txt 7.2 Botanical Simulator 76 4. Open the simulation configure file and modify the following lines (even better, copy the simulation configure file and edit the copy) : *domainList #base-domainList.txt base-myBigTree-domainList.forest.txt *parentList #base-parentList.txt base-myBigTree-parentList.forest.txt Here, # means that the line is disabled. In other words, we are disabling the old hierarchy and enabling a new one. The forest simulation can be launched at this point. Currently, only replicated copies of the same plant can be used. But this could be changed to allow different plants in the same scene. 7.2 Botanical Simulator Our real-time simulator must be launched with a simulation configure file generated by the preprocessor. This simulation file contains several initial settings, such as timestep, damping coefficients, number of threads, camera position and angles, window size, etc., and data file- names such as the domain list filename and parent list filename. Figure 7.6 gives a screenshot of the runtime driver. User is able to pull on the plant with the left mouse button, apply gravity (on GUI), or use Perlin wind (on GUI, or specified in configure file). 7.2.1 Simulation Control The simulation control panel is presented in Figure 7.7. Procedural damping User is able to control the damping via the dimensionlessξ parameter (see Section 5.3). 7.2 Botanical Simulator 77 Figure 7.6: Botanical Simulator: Left: Simulation display window. It allows user to interactive with the simulation by mouse. Middle: Frequency tuning panel. User could adjust the material stiffness on this panel. Right: Simulation control panel. User could change timestep, damping, mouse force strength, wind etc. on this panel. 1. ξ = 0: no damping 2. 0<ξ < 1: under-damped system 3. ξ = 1: critical damping 4. ξ > 1: over damped system User can control the value of ξ as a function of the frequency of each domain. This is done by modelingξ as a constant function ofν forν ≤ ν low andν ≥ ν high , and a straight line in between (as explained in Section 5.3). Quality Quality refers to the quality control coefficient introduced in Section 5.2. The sim- ulator uses asynchronous timestepping (see Section 5.4). It will step each domain at such a timestep that Q steps are made per longest oscillation period (lowest frequency) of this domain. Parameter Q (quality) therefore controls how many integration steps are performed, relative to the period of the domain. Q= 100 is a good initial value. Higher values of Q will produce 7.2 Botanical Simulator 78 better motion (less numerical damping), at the expense of more computation. The “Set to 1x speed” button applies a heuristic to tune quality in such a way that the simulation runs at the same rate as physics. Figure 7.7: Simulation control panel Timestep This sets the graphics timestep on a phys- ical scale. It is usually set to 1 / 30 s or 1 / 24 s. When Timestep quality control is enabled, asyn- chronous timestepping is used. The timestep for each domain is subdivided automatically, as explained in the “Quality” section above (also see Section 5.2). Liveliness and Inertia These two important param- eters control the strength of system forces (forces ap- plied to a branch because of the motion of its par- ent, Section 3.3) and inertia forces (forces applied to a branch because of the motion of its children, Section 3.3). Tuning these values well is important. Default liveliness is 1.0 but the user can increase liveliness to 2.0 or 3.0, causing the tree to be more lively. In such cases, Q (quality) may need to be increased somewhat. User can also decrease the liveliness to a value below 1.0. This may help stabilizing the simulation without increasing Q. Hierarchical levels: This makes it possible to sim- ulate only a subset of the tree (levels up to a given hierarchical depth). This is controlled via “Max simulated level”. You can also make levels beyond a certain level rigid, via “Max flexible level”. “Max simulated frequency” is a cutoff frequency, namely, domains above that frequency 7.2 Botanical Simulator 79 will be made rigid as they vibrate too quickly to be seen in the animation. Typically it is set to 20 Hz. Perlin wind See Chapter 6. User can enable it on the GUI. Standard Perlin wind parameters are used. User can specify a constant wind direction plus the space-time 4D Perlin wind. To visualize the Perlin wind, press ‘c’ to show the Perlin wind field, and/or ‘ C’ to show the Perlin wind forces applied on sample points. Animation player Makes it possible to record the simulation to a buffer and then play it back. Rendering We describe some rendering switches on keyboard here. Use ‘e’, ‘E’, ‘w’, ‘W’ to show or hide solid voxel mesh, solid triangle mesh, wireframe voxel mesh and wireframe triangle mesh. By pressing, ‘1’, ‘2’, ‘3’, user can show/hide branches, fronds and leaves. 7.2.2 Frequency Tuning This is an independent GUI panel to tune the aging mode and stiffness of each hierarchy level at runtime. See Figure 7.8. Frequency tuning This panel makes it possible to set the stiffness multiplicator individually for each individual level. A global multiplicator (multiplies all levels simultaneously, in ad- dition to the level multiplicators) is available via “Freq Scaling”. User can enable “aging mode” (procedural stiffness, Section 5.1) or “absolute mode” (constant stiffness). With pro- cedural stiffness, stiffness is determined based on age of each branch. Age is determined based on the branch geometry (size). Aging mode has the exponent parameter that controls how aggressive the aging mode is. A value of 2.0 is a good choice. Higher values will make small branches more lively. A value of 0.0 is equivalent to absolute mode (constant stiffness). Aging randomization randomizes the frequency somewhat, so that we don’t get an entire tree resonating at 1Hz in the aging mode with exponent 2. Frequencies will be multiplied with 7.2 Botanical Simulator 80 a random factor sampled uniformly on the interval [1.0/(1.0+ agingRandomization),1.0+ agingRandomization]. Note that user needs to select the levels before clicking on “Upload”. Important: aging mode or absolute mode changes will only apply to the selected levels. User can select all levels via “Select all”. A good way to tune a tree is to set the first few levels stiff, i.e., 3.0, 5.0 or 10.0. Then, make higher levels softer, i.e., 1.0. Figure 7.8: Frequency tuning panel Chapter 8 Conclusion We presented a robust system for stable physically based simulation of anatomically realistic botanical systems. We first presented a real-time algorithm for simulation of reduced nonlinear flexible multibody systems undergoing large deformations. The algorithm was made possible by deriving the first and second time derivatives of the rotation matrix used in polar decom- position. The algorithm supports localized deformations, requires no constraints, and runs in time linear in the number of domains. This algorithm forms the foundation of our botanical simulator. We then demonstrated a preprocessor to pre-process “polygon soup” plant geometry for do- main decomposition simulations. Our system scales to the complexity of real-world adult trees, flowers and bushes. We have pre-processed over 100 plants from several publicly avail- able vegetation model libraries. Our system accommodates unorganized, unprocessed triangle input geometry, including billboards. We extended our simulator to support fracture, interac- tive plant design, and forest simulation. Third, we proposed a series of procedural methods to automatically tune timestep, stiffness, and damping coefficient for every domain. These tuning techniques could be incorporated into any modal-based FEM simulation method. We extended our single-core multi-domain dy- namics algorithm to an asynchronous parallel timestepping algorithm that utilizes multi-core computing. We demonstrated that this algorithm greatly enhanced the simulation speed and numerical stability under large deformations. 82 Finally we described how to produce a robust, 4D randomized wind force field necessary for plant simulation. Limitations and future work: Our multi-domain dynamics algorithm is limited to domain topologies without loops. Several bodies connected and rooted to the ground form a loop, and must be simulated as a single domain in our system. Loops could be closed by adding springs; or more formally, using extensions paralleling those for rigid articulated chains [26]. We as- sume that the domain interfaces undergo only a small amount of non-rigid deformation. This assumption is valid when the interfaces are small, and worked well in our examples. Flexible interfaces could be simulated by adding additional “boundary” modes to each domain [92]. Related to that, it may seem plausible to add additional terms to the equations of motion that would explicitly couple the elastic deformations of two domains meeting at an interface. Note that for strictly rigid interfaces such terms are identically zero, and that any external forces al- ready are properly propagated to parent domains via our interface forces. In our work, frames follow the motion imposed by the parent exactly, without any freedom for deviation. Some of the six degrees of freedom could be relaxed by introducing joints, which would lead to a method supporting both articulation and large deformations. Our preprocessor cannot handle plants that are in continuous contact with their natural envi- ronment, such as a vine climbing a mesh fence or ivy climbing a tree. For plants in simpler (ground) contact such as zucchini or watermelons, frictional contact could be handled using constraint solvers, or even penalty forces. Our system avoids loops, which has not been a prob- lem in practice as the vast majority of plants do not have loops. Loops could be addressed using penalty forces. Our material parameters could in the future incorporate mechanical properties or other observation data from real plants. Our instancing can be easily extended to domains that consist of several disjoint components, each texture-mapped with a distinct texture map. A more challenging case, however, would be to support hierarchical instancing, where instances themselves consist of replicated copies, e.g., replicated blooms, each consisting of replicated but otherwise identical petals. We note, however, that many blooms in practice are modeled as a single global quadrilateral billboard, replicated in different orientations, where our single- 83 level instancing is sufficient. Methods that use alternative simulation bases [29] may be an interesting approach to simulate plants. We use a semi-implicit integrator to timestep our models which provides stability but also introduces artificial damping. It would be interesting to seek integration schemes that can work with model reduction and that can avoid artificial damping. Simulator realism would be improved by handling self-collisions, and simulating two-way coupling between the wind and the plant. Bibliography [1] Akagi, Y . and Kitajima, K. (2006). Computer animation of swaying trees based on physical simulation. Computers & Graphics, 30(4):529–539. [2] Almroth, B. O., Stern, P., and Brogan, F. A. (1978). Automatic Choice of Global Shape Functions in Structural Analysis. AIAA Journal, 16(5):525–528. [3] An, S. S., Kim, T., and James, D. L. (2008). Optimizing cubature for efficient integration of subspace deformations. ACM Trans. on Graphics, 27(5):165:1–165:10. [4] Anonymous (2012). Personal correspondence with a leading special effects studio. [5] Antoulas, A. C., Sorensen, D. C., and Gugercin, S. (2001). A survey of model reduction methods for large-scale systems. Contemporary Mathematics, 280:193–219. [6] Araldsen, P. O. (1974). An Example of Large-Scale Structural Analysis. Comparison between Finite Element Calculation and Full Scale Measurements on The Oil Tanker Esso Norway. Computers and Structures, 4(1):69–93. [7] Barbiˇ c, J., da Silva, M., and Popovi´ c, J. (2009). Deformable object animation using re- duced optimal control. ACM Trans. on Graphics, 28(3). [8] Barbiˇ c, J. and James, D. L. (2005). Real-time subspace integration for St. Venant- Kirchhoff deformable models. ACM Trans. on Graphics, 24(3):982–990. [9] Barbiˇ c, J. and Zhao, Y . (2011). Real-time large-deformation substructuring. ACM Trans. on Graphics (SIGGRAPH 2011), 30(4):91:1–91:7. [10] Beaudoin, J. and Keyser, J. (2004). Simulation levels of detail for plant motion. In Symp. on Computer Animation (SCA), pages 297–304. [11] Benner, P., Freund, R. W., Sorensen, D. C., and Varga, A. (2006). Special Issue on Order Reduction of Large-Scale Systems. Linear Algebra and its Applications, 415:231–234. [12] Bertails, F. (2009). Linear time super-helices. Comput. Graphics Forum, 28(2):417–426. [13] Bickel, B., Baecher, M., Otaduy, M., Matusik, W., Pfister, H., and Gross, M. (2009). Capture and modeling of non-linear heterogeneous soft tissue. ACM Trans. on Graphics (SIGGRAPH 2009), 28(3):89:1–89:9. Bibliography 85 [14] Bloomenthal, J. (1985). Modeling the mighty maple. In Computer Graphics (Proc. of SIGGRAPH 88), volume 19, pages 305–311. [15] Boyer, F. and Khalil, W. (1995). Recursive inverse and forward dynamics of flexible manipulators. In Proc. Conference ECC, Rome, Italy, pages 2696–2701. [16] Capell, S., Green, S., Curless, B., Duchamp, T., and Popovi´ c, Z. (2002a). A Multireso- lution Framework for Dynamic Deformations. In Proc. of the Symp. on Comp. Animation 2002, pages 41–48. [17] Capell, S., Green, S., Curless, B., Duchamp, T., and Popovi´ c, Z. (2002b). Interactive skeleton-driven dynamic deformations. ACM Trans. on Graphics, 21(3):586–593. [18] Changizi, K. and Shabana, A. A. (1988). A recursive formulation for the dynamic analy- sis of open loop deformable multibody systems. Journal of Applied Mechanics, 55(3):687– 693. [19] Chao, I., Pinkall, U., Sanan, P., and Schröder, P. (2010). A Simple Geometric Model for Elastic Deformations. ACM Transactions on Graphics, 29(3):38:1–38:6. [20] Cormen, T. H., Leiserson, C. E., and Rivest, R. L. (1990). Introduction to Algorithms. MIT Press/McGraw-Hill. [21] Craig, R. and Bampton, M. (1968). Coupling of substructures for dynamic analysis. AIAA Journal, 6(7):1313–1319. [22] Debunne, G., Desbrun, M., Cani, M.-P., and Barr, A. H. (2001). Dynamic Real-Time Deformations Using Space & Time Adaptive Sampling. In Proc. of ACM SIGGRAPH 2001, pages 31–36. [23] Deussen, O. and Lintermann, B. (2005). Digital Design of Nature: Computer Generated Plants and Organics. Springer-Verlag, New York. [24] Diener, J., Rodriguez, M., Baboud, L., and Reveret, L. (2009). Wind projection basis for real-time animation of trees. Computer Graphics Forum, 28(2):533–540. [25] Dodds, R. H. and Lopez, L. A. (1980). Substructuring in Linear and Nonlinear Analysis. International Journal for Numerical Methods in Engineering, 15:583–597. [26] Featherstone, R. (1987). Robot Dynamics Algorithms. Kluwer Academic Publishers, Boston. [27] Georgii, J. and Westermann, R. (2005). A multigrid framework for real-time simulation of deformable volumes. In Proc. of the 2nd Workshop On Virtual Reality Interaction and Physical Simulation, pages 50–57. [28] Gildin, E. (2006). Model and controller reduction of large-scale structures based on projection methods. PhD thesis, The Institute for Computational Engineering and Sciences, University of Texas at Austin. Bibliography 86 [29] Gilles, B., Bousquet, G., Faure, F., and Pai, D. K. (2011). Frame-based elastic models. ACM Trans. on Graphics, 30(2):15:1–15:12. [30] Grimme, E., Sorensen, D., and Dooren, P. V . (1995). Model reduction of state space systems via an implicitly restarted Lanczos method. Numerical Algorithms, 12:1–31. [31] Grinspun, E., Krysl, P., and Schröder, P. (2002). CHARMS: A Simple Framework for Adaptive Simulation. In Proc. of ACM SIGGRAPH 2002. [32] Habel, R., Kusternig, A., and Wimmer, M. (2009). Physically Guided Animation of Trees. Computer Graphics Forum, 28(2):523–532. [33] Hay, A., Borggaard, J. T., and Pelletier, D. (2009). Local Improvements to Reduced- Order Models Using Sensitivity Analysis of the Proper Orthogonal Decomposition. J. of Fluid Mechanics, 629:41–72. [34] Hu, S., Chiba, N., and He, D. (2012). Realistic animation of interactive trees. The Visual Computer, 28:859–868. [35] Huang, J., Liu, X., Bao, H., Guo, B., and Shum, H.-Y . (2006). An efficient large defor- mation method using domain decomposition. Computers & Graphics, 30(6):927 – 935. [36] Ijiri, T., Owada, S., and Igarashi, T. (2006). The sketch L- System: Global control of tree modeling using free-form strokes. In Smart Graphics, pages 138–146. [37] Interactive Data Visualization (1999). Speedtree. www.speedtree.com. [38] J. V . Clark and N. Zhou and D. Bindel and L. Schenato and W. Wu and J. Demmel and K. S. J. Pister (2000). 3d mems simulation using modified nodal analysis. In Proc. of Microscale Systems: Mechanics and Measurements Symposium, pages 68–75. [39] James, D. L., Barbiˇ c, J., and Twigg, C. D. (2004). Squashing Cubes: Automating De- formable Model Construction for Graphics. In Proc. of ACM SIGGRAPH Sketches and Applications. [40] James, D. L. and Pai, D. K. (2002a). DyRT: Dynamic Response Textures for Real Time Deformation Simulation With Graphics Hardware. ACM Trans. on Graphics, 21(3):582– 585. [41] James, D. L. and Pai, D. K. (2002b). Real Time Simulation of Multizone Elastokinematic Models. In IEEE Int. Conf. on Robotics and Automation, pages 927–932. [42] James, D. L., Twigg, C. D., Cove, A., and Wang, R. Y . (2007). Mesh Ensemble Motion Graphs: Data-driven Mesh Animation with Constraints. ACM Trans. on Graphics, 26(4). [43] James, K. R., Haritos, N., and Ades, P. K. (2006). Mechanical stability of trees under dynamic loads. American J. of Botany, 93(10):1522–1530. Bibliography 87 [44] Kaufman, D. M., Sueda, S., James, D. L., and Pai, D. K. (2008). Staggered Projections for Frictional Contact in Multibody Systems. ACM Transactions on Graphics, 27(5):164:1– 164:11. [45] Kaufmann, P., Martin, S., Botsch, M., and Gross, M. (2009). Flexible Simulation of De- formable Models Using Discontinuous Galerkin FEM. J. of Graphical Models, 71(4):153– 167. [46] Kharevych, L., Mullen, P., Owhadi, H., and Desbrun, M. (2009). Numerical coarsening of inhomogeneous elastic materials. ACM Trans. on Graphics, 28(3):51:1–51:8. [47] Kim, T. and James, D. (2009a). Skipping steps in deformable simulation with online model reduction. ACM Trans. on Graphics (SIGGRAPH Asia 2009), 28(5):123:1–123:9. [48] Kim, T. and James, D. (2011). Physics-based character skinning using multi-domain subspace deformations. In Symp. on Computer Animation (SCA), pages 63–72. [49] Kim, T. and James, D. L. (2009b). Skipping steps in deformable simulation with online model reduction. ACM Trans. Graph., 28(5):1–9. [50] Ko, J., Kurdila, A. J., and Rediniotis, O. K. (2000). Divergence free bases and multires- olution methods for reduced-order flow modeling. AIAA Journal, 38(12):2219–2232. [51] Krysl, P., Lall, S., and Marsden, J. E. (2001). Dimensional model reduction in non- linear finite element dynamics of solids and structures. Int. J. for Numerical Methods in Engineering, 51:479–504. [52] Labelle, F. and Shewchuk, J. R. (2007). Isosurface Stuffing: Fast Tetrahedral Meshes with Good Dihedral Angles. In Proc. of ACM SIGGRAPH 2007, pages 471–478. [53] Lall, S., Marsden, J. E., and Glavaski, S. (2002). A subspace approach to balanced truncation for model reduction of nonlinear control systems. Int. J. on Robust and Nonlinear Control, 12:519–535. [54] Lehoucq, R., Sorensen, D., and Yang, C. (1997). ARPACK Users’ Guide: Solution of large scale eigenvalue problems with implicitly restarted Arnoldi methods. Technical report, Comp. and Applied Mathematics, Rice Univ. [55] Li, R.-C. and Bai, Z. (2005). Structure preserving model reduction using a Krylov sub- space projection formulation. Comm. Math. Sci., 3(2):179–199. [56] Lin, M. C. and Gottschalk, S. (1998). Collision Detection Between Geometric Models: A Survey. In Proc. of IMA Conference on Mathematics of Surfaces, pages 37–56. [57] Lindenmayer, A. (1968). Mathematical models for cellular interaction in development. J. of Theoretical Biology Parts I and II, 18:280–315. [58] Lintermann, B. and Deussen, O. (1999). Interactive modeling of plants. IEEE Computer Graphics and Applications, 19(1):56–65. Bibliography 88 [59] Lu, H., Guo, X., Zhao, C., and Li, C. (2011). Physical model for interactive deformation of 3d plant. International Journal of Virtual Reality, 10(2):33. [60] Lumley, J. L. (1967). The structure of inhomogeneous turbulence. In A.M.Yaglom and V .I.Tatarski, editors, Atmospheric turbulence and wave propagation, pages 166–178. [61] Mathai, A. M. (1997). Jacobians of Matrix Transformations and Functions of Matrix Argument. World Scientific Publishing Co. [62] Metaxas, D. and Terzopoulos, D. (1992). Dynamic deformation of solid primitives with constraints. In Computer Graphics (Proc. of SIGGRAPH 92), pages 309–312. [63] Mezger, J., Thomaszewski, B., Pabst, S., and Strasser, W. (2008). Interactive physically- based shape editing. In Proc. of the ACM Symp. on Solid and physical modeling, pages 79–89. [64] Moore, B. C. (1981). Principal Component Analysis in Linear System: Controllability, Observability and Model Reduction. IEEE Trans. on Automatic Control, 26:17–32. [65] Müller, M. and Gross, M. (2004). Interactive Virtual Materials. In Proc. of Graphics Interface 2004, pages 239–246. [66] Müller, M., Heidelberger, B., Teschner, M., and Gross, M. (2005). Meshless Deforma- tions Based on Shape Matching. In Proc. of ACM SIGGRAPH 2005, pages 471–478. [67] Mˇ ech, R. and Prusinkiewicz, P. (1996). Visual models of plants interacting with their environment. In Proc. of ACM SIGGRAPH 1996, pages 397–410. [68] Nealen, A., Müller, M., Keiser, R., Boxerman, E., and Carlson, M. (2006). Physically based deformable models in computer graphics. Computer Graphics Forum, 25(4):809– 836. [69] Nesme, M., Kry, P. G., Jeˇ rábková, L., and Faure, F. (2009). Preserving topology and elasticity for embedded deformable models. ACM Trans. on Graphics, 28(3):52:1–52:9. [70] Nickell, R. E. (1976). Nonlinear Dynamics by Mode Superposition. Computer Methods in Applied Mechanics and Engineering, 7:107–129. [71] Noor, A. K., Kamel, H. A., and Fulton, R. E. (1977). Substructuring Techniques - Status and Projections. Computer and structures, 8:621–632. [72] Oppenheimer, P. E. (1986). Real time design and animation of fractal plants and trees. In Proc. of ACM SIGGRAPH 1986, pages 55–64. [73] Ota, S., Tamura, M., Fujimoto, T., Muraoka, K., and Chiba, N. (2004). A hybrid method for real-time animation of trees swaying in wind fields. The Visual Computer, 20:613–623. [74] Parker, E. G. and O’Brien, J. F. (2009). Real-time deformation and fracture in a game environment. In Symp. on Computer Animation (SCA), pages 156–166. Bibliography 89 [75] Patnaik, S., Gendy, A., and Hopkins, D. (1994). Design optimization of large structural systems with substructuring in a parallel computational environment. Computing Systems in Engineering, 5(4-6):425–440. [76] Pentland, A. and Williams, J. (1989). Good vibrations: Modal dynamics for graphics and animation. Computer Graphics (Proc. of ACM SIGGRAPH 89), 23(3):215–222. [77] Perlin, K. (2002). Improving Noise. ACM Trans. on Graphics, 21(3):681–682. [78] Pirk, S., Stava, O., Kratt, J., Said, M. A. M., Neubert, B., Mˇ ech, R., Benes, B., and Deussen, O. (2012). Plastic trees: interactive self-adapting botanical tree models. ACM Trans. Graph., 31(4):50:1–50:10. [79] Prusinkiewicz, P. (1986). Graphical applications of l-systems. In Graphics Interface / Vision Interface, pages 247–253. [80] Reche-Martinez, A., Martin, I., and Drettakis, G. (2004). V olumetric reconstruction and interactive rendering of trees from photographs. In Proc. of ACM SIGGRAPH 2004, pages 720–727. [81] Rewienski, M. and White, J. (2006). Model order reduction for nonlinear dynamical systems based on trajectory piecewise-linear approximations. Linear Algebra and its Ap- plications, 412(2-3):426–454. [82] Ryu, Y . S. and Arora, J. S. (1985). Review of Nonlinear FE Methods with Substructures. Journal of Engineering Mechanics, 111(11):1361–1379. [83] Sakaguchi, T. and Ohya, J. (1999). Modeling and animation of botanical trees for interac- tive virtual environments. In Proc. of the Symp. on Virtual reality software and technology, pages 139–146. [84] Sellgren, U. (2003). Component Mode Synthesis - A Method for Efficient Dynamic Sim- ulation of Complex Technical Systems. Technical report, The Royal Institute of Technology (KTH), Stockholm, Sweden. [85] Shabana, A. A. (2005). Dynamics of Multibody Systems. Cambridge Univ. Press, New York, NY. [86] Sharf, I. and D’Eleuterio, G. (1988). Computer simulation of elastic chains using a recursive formulation. In IEEE Conf. on Robotics and Automation, volume 3, pages 1539– 1545. [87] Sifakis, E. and Barbiˇ c, J. (2012). FEM Simulation of 3D Deformable Solids: A prac- titioner’s guide to theory, discretization and model reduction, Part 2: Model reduction. In SIGGRAPH Course Notes. www.femdefo.org. [88] Simeoni, M., Vandenbosch, G. A. E., and Lager, I. E. (2005). Model-order reduction techniques for linear electromagnetic problems - an overview. In European Microwave Conference, page 4 pp. Bibliography 90 [89] Simon Green (2005). Implementing Improved Perlin Noise. In Matt Pharr, editor, GPU Gems 2, pages 409–416. Addison-Wesley. [90] Sousa, T. and Crytek (2007). GPU Gems 3, Chapter 16. Vegetation Procedural Animation and Shading in Crysis. Addison-Wesley Professional, Boston. [91] Stam, J. (1997). Stochastic Dynamics: Simulating the Effects of Turbulence on Flexible Structures. Comp. Graphics Forum, 16(3). [92] Storaasli, O. O. and Bergan, P. (1987). Nonlinear Substructuring Method for Concurrent Processing Computers. AIAA, 25(6):871–876. [93] Toselli, A. and Widlund, O. (2004). Domain Decomposition Methods. Springer. [94] Treuille, A., Lewis, A., and Popovi´ c, Z. (2006). Model reduction for real-time fluids. ACM Trans. on Graphics, 25(3):826–834. [95] Turbosquid (2000). www.turbosquid.com. [96] Twigg, C. and Kaˇ ci´ c-Alesi´ c, Z. (2010). Point cloud glue: constraining simulations using the procrustes transform. In Symp. on Computer Animation (SCA), pages 45–54. [97] Twigg, C. and Kaˇ ci´ c-Alesi´ c, Z. (2011). Optimization for sag-free simulations. In Symp. on Computer Animation (SCA), pages 225–236. [98] Weber, J. P. (2008). Fast simulation of realistic trees. IEEE Computer Graphics and Applications, 28(3):67–75. [99] Wicke, M., Stanton, M., and Treuille, A. (2009). Modular bases for fluid dynamics. ACM Trans. on Graphics, 28(3):39:1–39:8. [100] Wong, J. C. and Datta, A. (2004). Animating real-time realistic movements in small plants. In Proc. of GRAPHITE 2004, pages 182–189. [101] Xfrog (2009). www.xfrog.com. [102] Yang, Y ., Xu, W., Guo, X., Zhou, K., and Guo, B. (2012). Boundary-aware multi- domain subspace deformation. IEEE Trans. on Visualization and Computer Graphics, pending minor revision. http://www.utdallas.edu/˜yxy061100/. [103] Zhang, L., Song, C., Tan, Q., Chen, W., and Peng, Q. (2006). Quasi-physical Simulation of Large-Scale Dynamic Forest Scenes, volume 4035 of Lecture Notes in Computer Science. Springer. [104] Zhang, L., Zhang, Y ., Chen, W., and Peng, Q. (2008). Real-time simulation of large- scale dynamic forest with gpu. In IEEE Conf. on Circuits and Systems, page 614–617. [105] Zhang, L., Zhang, Y ., Jiang, Z., Li, L., Chen, W., and Peng, Q. (2007). Precomputing data-driven tree animation. Computer Animation and Virtual Worlds, 18(4-5):371–382. [106] Zhao, Y . and Barbiˇ c, J. (2013). Interactive authoring of simulation-ready plants. ACM Trans. on Graphics (SIGGRAPH 2013), 32(4):84:1–84:12.

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Anatomically based human hand modeling and simulation

PDF

CUDA deformers for model reduction

PDF

Acquisition of human tissue elasticity properties using pressure sensors

PDF

Artistic control combined with contact and elasticity modeling in computer animation pipelines

PDF

Closing the reality gap via simulation-based inference and control

PDF

Multi-scale dynamic capture for high quality digital humans

PDF

Rethinking perception-action loops via interactive perception and learned representations

PDF

High-throughput methods for simulation and deep reinforcement learning

PDF

Adaptive event-driven simulation strategies for accurate and high performance retinal simulation

PDF

3D deep learning for perception and modeling

PDF

Feature-preserving simplification and sketch-based creation of 3D models

PDF

Deformation control for mask image projection based stereolithography process

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Point-based representations for 3D perception and reconstruction

PDF

Data-driven acquisition of closed-loop robotic skills

PDF

Dynamic modeling and simulation of flapping-wing micro air vehicles

PDF

Modeling and simulation of complex recovery processes

PDF

Modeling information operations and diffusion on social media networks

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Model based design of porous and patterned surfaces for passive turbulence control

Asset Metadata

Creator Zhao, Yili (author)

Core Title Plant substructuring and real-time simulation using model reduction

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Computer Science

Publication Date 08/20/2014

Defense Date 04/17/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag authoring,botanical,computer animation,computer graphics,domain decomposition,FEM,interactive,large deformations,model reduction,nonlinear elasticity,OAI-PMH Harvest,physically‐based simulation,plant,substructuring

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Barbič, Jernej (committee chair), Kukavica, Igor (committee member), Neumann, Ulrich (committee member), Schaal, Stefan (committee member), Sukhatme, Gaurav S. (committee member)

Creator Email yilizhao.pku@gmail.com,yilizhao@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-462785

Unique identifier UC11287579

Identifier etd-ZhaoYili-2834.pdf (filename),usctheses-c3-462785 (legacy record id)

Legacy Identifier etd-ZhaoYili-2834.pdf

Dmrecord 462785

Document Type Dissertation

Format application/pdf (imt)

Rights Zhao, Yili

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA