Science, Faculty of
Computer Science, Department of
DSpace
UBCV
Zhu, Yufeng
Doctor of Philosophy - PhD
University of British Columbia
This thesis investigates several aspects of the optimal mapping problem, finding a bijective function between two shapes which minimizes some metric, and its many applications in computer graphics such as finding planar embeddings of curved surface meshes, image warping, mesh deformation, elastoplastic simulation, and more. We highlight in particular problems where differences in topology between the shapes—which necessitate discontinuities in the mapping—play a crucial role, e.g. due to animation decisions or physical fracture. Our first project presents an approach to extend discrete variational shape interpolation to create keyframe animation involving extreme deformation, topology change and dynamics. To construct correspondence between keyframe shapes as well as satisfying feature matching constraints, we introduce a consistent multimesh structure that is able to resolve topological in-equivalence between shapes and a Comesh Optimization algorithm that optimizes our multimesh for both intra and inter-mesh quality. To further speed up the total solving time, we then develop three new improvements to the state of the art nonlinear optimization techniques that are frequently used in this field, including a barrier-aware line-search filter that cures blocked descent steps due to element barrier terms; an energy proxy model that adaptively blends the Sobolev gradient and L-BFGS descent to gain the advantages of both; and a characteristic gradient norm providing a robust and largely mesh- and energy independent convergence criterion. Finally, we demonstrate another related interesting graphics application, brittle fracture simulation. In particular, we investigate whether we can sidestep the volumetric optimal mapping problem for solving elastic deformation and instead use a boundary element discretization which only requires a surface mesh. We will show how the computational cost of such problems can be alleviated using a boundary integral formulation and kernel-independent Fast Multipole Method. By combining with an explicit surface tracking framework, we further avoid the expensive volumetric mesh construction and maintenance during fracture propagation and shape topology change.
Optimal Mapping with Topology Change for Animationand Geometry ProblemsbyYufeng ZhuB. Sc, Shanghai Jiaotong University, 2010M. Sc, University of Southern California, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University of British Columbia(Vancouver)December 2018c© Yufeng Zhu, 2018The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Optimal Mapping with Topology Change for Animation and GeometryProblemssubmitted by Yufeng Zhu in partial fulfillment of the requirements for the degreeof Doctor of Philosophy in Computer Science.Examining Committee:Robert Bridson, Computer ScienceSupervisorChen Greif, Computer ScienceCo-supervisorAlla Sheffer, Computer ScienceUniversity ExaminerCarl Ollivier-Gooch, Mechanical EngineeringUniversity ExaminerPaul Kry, Computer Science, McGill UniversityExternal ExamineriiAbstractThis thesis investigates several aspects of the optimal mapping problem, finding abijective function between two shapes which minimizes some metric, and its manyapplications in computer graphics such as finding planar embeddings of curvedsurface meshes, image warping, mesh deformation, elastoplastic simulation, andmore. We highlight in particular problems where differences in topology betweenthe shapes—which necessitate discontinuities in the mapping—play a crucial role,e.g. due to animation decisions or physical fracture.Our first project presents an approach to extend discrete variational shape in-terpolation to create keyframe animation involving extreme deformation, topologychange and dynamics. To construct correspondence between keyframe shapes aswell as satisfying feature matching constraints, we introduce a consistent multi-mesh structure that is able to resolve topological in-equivalence between shapesand a Comesh Optimization algorithm that optimizes our multimesh for both intra-and inter-mesh quality.To further speed up the total solving time, we then develop three new improve-ments to the state of the art nonlinear optimization techniques that are frequentlyused in this field, including a barrier-aware line-search filter that cures blocked de-scent steps due to element barrier terms; an energy proxy model that adaptivelyblends the Sobolev gradient and L-BFGS descent to gain the advantages of both;and a characteristic gradient norm providing a robust and largely mesh- and energy-independent convergence criterion.Finally, we demonstrate another related interesting graphics application, brittlefracture simulation. In particular, we investigate whether we can sidestep the volu-metric optimal mapping problem for solving elastic deformation and instead use aiiiboundary element discretization which only requires a surface mesh. We will showhow the computational cost of such problems can be alleviated using a boundaryintegral formulation and kernel-independent Fast Multipole Method. By combin-ing with an explicit surface tracking framework, we further avoid the expensivevolumetric mesh construction and maintenance during fracture propagation andshape topology change.ivLay SummaryWe may learn about set and mapping in our high school’s math class. In short, aset is to put things or elements with similar properties into an abstract category.For example, among a group of people, men and women can be grouped into twosets. Mapping is the way to describe the relation between sets of elements. In theprevious example, mapping can be defined as couple relationship. Such conceptsmay sounds abstract and boring but they are widely used in our daily life, especiallyin scientific research and industrial application. In this work, I studied ways to findoptimal structure preserving map with their applications in computer graphics field,such as texture mapping, 3D geometry modeling, elastoplastic simulation, shapeinterpolation, etc. Such applications play crucial roles in industry areas, like visualeffect, computer games, virtual reality and so on.vPrefaceThe work presented in this thesis is co-supervised by Dr. Robert Bridson and Dr.Chen Greif at the University of British Columbia. It incorporates three publishedpapers as separate chapters:• Chapter 2 was published as Y. Zhu, J. Popovic´, R. Bridson and D. Kaufman.Planar Interpolation with Extreme Deformation, Topology Change and Dy-namics[157]. ACM Transactions on Graphics, 36(6), 2017 and presented atSIGGRAPH Asia 2017.• Chapter 3 was published as Y. Zhu, R. Bridson and D. Kaufman. BlendedCured Quasi-Newton for Distortion Optimization[158]. ACM Transactionson Graphics, 37(4), 2018 and presented at SIGGRAPH 2018.• Chapter 4 was published as Y. Zhu, R. Bridson and C. Greif. Simulat-ing Rigid Body Fracture with Surface Meshes[156]. ACM Transactions onGraphics, 34(4), 2015 and presented at SIGGRAPH 2015.Chapter 2 and Chapter 3 come from two summer research internship projects in theyears of 2015 and 2016 respectively at Creative Technology Lab, Adobe Seattle inUnited States. Dr. Danny Kaufman and Dr. Jovan Popovic´ were my internshipmentors in 2015 and Dr. Kaufman also mentored my internship project in 2016.Research directions were suggested by my supervisors, while the original researchideas were proposed by me and developed in concert with my collaborators.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Planar Interpolation with Extreme Deformation, Topology Changeand Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2.1 Planar Animation . . . . . . . . . . . . . . . . . . . . . . 72.2.2 Shape Interpolation . . . . . . . . . . . . . . . . . . . . . 72.2.3 Compatible Meshing . . . . . . . . . . . . . . . . . . . . 92.3 Planar Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . 92.3.1 Inbetweening via Variational Interpolation . . . . . . . . . 102.3.2 Multimesh . . . . . . . . . . . . . . . . . . . . . . . . . 12vii2.4 Comesh Optimization . . . . . . . . . . . . . . . . . . . . . . . . 152.4.1 Mesh and Mapping Energies . . . . . . . . . . . . . . . . 172.4.2 Optimizing Solely Mesh or Mapping is Insufficient . . . . 182.4.3 Multishape MIPS and Symmetry . . . . . . . . . . . . . . 192.4.4 Comeshing Constraints . . . . . . . . . . . . . . . . . . . 212.4.5 Comesh Optimization . . . . . . . . . . . . . . . . . . . 222.4.6 Solving Comesh Optimization . . . . . . . . . . . . . . . 232.5 Interpolation Energies . . . . . . . . . . . . . . . . . . . . . . . . 262.5.1 Cohesive Energy . . . . . . . . . . . . . . . . . . . . . . 272.5.2 Non-overlap and Deformation Energies . . . . . . . . . . 302.6 Animating Inbetweens . . . . . . . . . . . . . . . . . . . . . . . 312.6.1 Dynamic Effects for Variational Interpolation . . . . . . . 312.7 Experiments and Evaluations . . . . . . . . . . . . . . . . . . . . 332.7.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . 332.7.2 Comeshing . . . . . . . . . . . . . . . . . . . . . . . . . 342.7.3 Animation . . . . . . . . . . . . . . . . . . . . . . . . . 353 Blended Cured Quasi-Newton for Distortion Optimization . . . . . 433.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.2 Problem Statement and Overview . . . . . . . . . . . . . . . . . 453.2.1 Iterative Solvers for Nonlinear Minimization . . . . . . . 453.3 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3.1 Energies and Applications . . . . . . . . . . . . . . . . . 483.3.2 Energy Approximations . . . . . . . . . . . . . . . . . . 493.3.3 Line Search . . . . . . . . . . . . . . . . . . . . . . . . . 523.3.4 Termination . . . . . . . . . . . . . . . . . . . . . . . . . 533.4 Blended Quasi-Newton . . . . . . . . . . . . . . . . . . . . . . . 543.4.1 Greedy Laplacian Blending . . . . . . . . . . . . . . . . 563.4.2 Blended Quasi-Newton . . . . . . . . . . . . . . . . . . . 583.5 Barrier-Aware Line Search Filtering . . . . . . . . . . . . . . . . 583.5.1 Curing Line Search . . . . . . . . . . . . . . . . . . . . . 593.5.2 One-sided Barriers in Distortion Optimization . . . . . . . 603.5.3 Iterating Away from Collapse . . . . . . . . . . . . . . . 60viii3.5.4 Line Search Filtering . . . . . . . . . . . . . . . . . . . . 613.6 Termination Criteria . . . . . . . . . . . . . . . . . . . . . . . . . 623.6.1 Characteristic Gradient Norm . . . . . . . . . . . . . . . 633.7 The BCQN Algorthim . . . . . . . . . . . . . . . . . . . . . . . 653.7.1 Comparison with Other Algorithms . . . . . . . . . . . . 653.8 Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 673.8.1 Implementation . . . . . . . . . . . . . . . . . . . . . . . 673.8.2 Termination . . . . . . . . . . . . . . . . . . . . . . . . . 693.8.3 Newton-type Methods . . . . . . . . . . . . . . . . . . . 713.8.4 A Note on Solving Proxies and Pardiso . . . . . . . . . . 733.8.5 First-order Methods . . . . . . . . . . . . . . . . . . . . 733.8.6 Across-the-board Comparisons . . . . . . . . . . . . . . . 773.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 773.9.1 Limitations and Future Work . . . . . . . . . . . . . . . . 804 Simulating Rigid Body Fracture with Surface Meshes . . . . . . . . 824.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.3 Algorithm Overview . . . . . . . . . . . . . . . . . . . . . . . . 884.4 Boundary Element Simulation . . . . . . . . . . . . . . . . . . . 904.4.1 Elastostatic Constitutive Model . . . . . . . . . . . . . . 904.4.2 Indirect Boundary Integral . . . . . . . . . . . . . . . . . 914.4.3 Boundary Element Discretization . . . . . . . . . . . . . 934.5 Fast Multipole Method . . . . . . . . . . . . . . . . . . . . . . . 964.6 Fracture Propagation Tracking . . . . . . . . . . . . . . . . . . . 994.6.1 Mesh-Based Surface Tracking . . . . . . . . . . . . . . . 1004.6.2 Contact Force Model and Fracture Criteria . . . . . . . . 1014.6.3 Fracture Initiation and Propagation . . . . . . . . . . . . 1014.6.4 Remeshing Operation . . . . . . . . . . . . . . . . . . . . 1054.7 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . 1064.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109ix5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1105.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1105.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114A Planar Interpolation with Extreme Deformation, Topology Changeand Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129A.1 Constructing Multimesh Structure . . . . . . . . . . . . . . . . . 129A.2 Directed Graph Flow Implementation . . . . . . . . . . . . . . . 132A.3 Comesh Optimization Energy Gradient . . . . . . . . . . . . . . . 133A.4 Continuous Energies for Inbetweening . . . . . . . . . . . . . . . 135A.5 Cohesive Energy Discretization . . . . . . . . . . . . . . . . . . . 137A.6 Transitivity of MIPS . . . . . . . . . . . . . . . . . . . . . . . . 139A.7 Automated Blending of Dynamic Effects . . . . . . . . . . . . . . 142A.8 Animation Results . . . . . . . . . . . . . . . . . . . . . . . . . 143A.9 Experiment Statistics . . . . . . . . . . . . . . . . . . . . . . . . 145A.9.1 Relative Weighting of Energies in Comesh OptimizationObjective . . . . . . . . . . . . . . . . . . . . . . . . . . 145A.9.2 Mesh Resolution . . . . . . . . . . . . . . . . . . . . . . 146A.9.3 Number of Example Shapes . . . . . . . . . . . . . . . . 147A.9.4 Problem Difficulty . . . . . . . . . . . . . . . . . . . . . 148B Blended Cured Quasi-Newton for Distortion Optimization . . . . . 149B.1 Equivalence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149C Simulating Rigid Body Fracture with Surface Meshes . . . . . . . . 152C.1 Boundary Element Method . . . . . . . . . . . . . . . . . . . . . 152C.1.1 Indirect Boundary Integral Formulation . . . . . . . . . . 152C.1.2 Potential Theory . . . . . . . . . . . . . . . . . . . . . . 152C.1.3 Elastostatic Model . . . . . . . . . . . . . . . . . . . . . 155C.1.4 Discretization . . . . . . . . . . . . . . . . . . . . . . . . 156C.2 Fast Multipole Method . . . . . . . . . . . . . . . . . . . . . . . 159C.2.1 Algorithm Introduction . . . . . . . . . . . . . . . . . . . 160xC.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166C.3.1 On BEM . . . . . . . . . . . . . . . . . . . . . . . . . . 166C.3.2 On FMM . . . . . . . . . . . . . . . . . . . . . . . . . . 168C.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . 172xiList of TablesTable 2.1 Math notation used for planar animation project . . . . . . . . 15Table 2.2 Comesh Optimization statistics of accompanying video examples 39Table 2.3 Animation statistics of accompany video examples . . . . . . . 40Table 4.1 Elastic deformation solved by the indirect boundary integralformulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Table 4.2 Detailed input information, algorithm timing and accuracy forscaling experiment . . . . . . . . . . . . . . . . . . . . . . . . 98Table 4.3 Runtime cost details for the 3D fracture simulation . . . . . . . 108Table A.1 Evaluation of relative contribution from the ODT mesh and cyclicMIPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145Table A.2 Evaluation of the Comesh Optimization’s scaling behavior bymesh refinement . . . . . . . . . . . . . . . . . . . . . . . . . 146Table A.3 Evaluation of the Comesh Optimization’s scaling behavior byincreasing meshes . . . . . . . . . . . . . . . . . . . . . . . . 147Table A.4 Evaluation of the Comesh Optimization’s performance by in-creasing the problem complexity . . . . . . . . . . . . . . . . 148Table C.1 Elastic deformation deformation solved by indirect boundaryintegral formulation . . . . . . . . . . . . . . . . . . . . . . . 174Table C.2 More elastic deformation solved by indirect boundary integralformulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 175xiiList of FiguresFigure 1.1 A UV parameterization example . . . . . . . . . . . . . . . . 2Figure 1.2 Topologically inequivalent shapes have mapping singularityproblem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2Figure 2.1 Keyframe interpolation examples with extreme deformationand topology change . . . . . . . . . . . . . . . . . . . . . . 5Figure 2.2 Square to arch example with different deformation path . . . . 8Figure 2.3 Our keyframe interpolation pipeline . . . . . . . . . . . . . . 10Figure 2.4 Unoptimized and Comesh Optimized keyframe interpolationcomparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 12Figure 2.5 Two topological operations to convert shapes into the sameequivalence class . . . . . . . . . . . . . . . . . . . . . . . . 13Figure 2.6 Cohesive edge pair illustration example . . . . . . . . . . . . 14Figure 2.7 A didactic example showing the benefit of Comesh Optimization 16Figure 2.8 Failure example of adopting previous meshing approaches . . 17Figure 2.9 Illustration of null space in piece wise linear discretization ofconformal mapping . . . . . . . . . . . . . . . . . . . . . . . 19Figure 2.10 Comesh Optimization result of a cell example . . . . . . . . . 21Figure 2.11 Illustration of cohesive edge constraints . . . . . . . . . . . . 22Figure 2.12 Illustration of topology update operations used in Comesh Op-timization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 2.13 Optimization result comparison of a square hexagon example 25Figure 2.14 Square to polygon interpolation example with optional sepa-rating control and dynamics . . . . . . . . . . . . . . . . . . 27xiiiFigure 2.15 Illustration of cohesive boundary potential kernel . . . . . . . 29Figure 2.16 Circle in circle example showing the ability of non-overlappinganimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31Figure 2.17 Timing statistics of Comesh Optimization scalability . . . . . 34Figure 2.18 Comesh Optimization evaluation experiment with incresinglyhard problems . . . . . . . . . . . . . . . . . . . . . . . . . . 36Figure 2.19 Trumperfly example showing the ability of our animation system 37Figure 2.20 Evaluation examples of Comesh Optimization algorithm’s scal-ability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38Figure 3.1 Optimization solver performance comparison over a twistingbar problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Figure 3.2 Illustration of severe effect on reducing step size from poorsearch direction . . . . . . . . . . . . . . . . . . . . . . . . . 47Figure 3.3 Illustration of problem dependency caused by standard termi-nation criterion . . . . . . . . . . . . . . . . . . . . . . . . . 53Figure 3.4 Comparison between unreliable stopping criterion and ours . . 54Figure 3.5 Performance comparison among various solvers on a planarconformal deformation problem . . . . . . . . . . . . . . . . 56Figure 3.6 Comparison among search direction curing strategies . . . . . 59Figure 3.7 Performance comparison among various solvers on a UV pa-rameterization problem . . . . . . . . . . . . . . . . . . . . . 63Figure 3.8 Sparsity pattern comparison between Laplacian proxy and otherproxies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 3.9 Convergence comparison of different termination criteria . . . 69Figure 3.10 Measuring improvement comparison of different terminationcriteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70Figure 3.11 Termination criteria comparison across mesh refinement andscale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71Figure 3.12 UV parameterization scaling, timing and sparsity comparisonamong different solvers . . . . . . . . . . . . . . . . . . . . . 72Figure 3.13 3D deformation scaling, timing and sparsity comparison amongdifferent solvers . . . . . . . . . . . . . . . . . . . . . . . . . 74xivFigure 3.14 Performance comparison between BCQN and PN using ar-madillo deformation test . . . . . . . . . . . . . . . . . . . . 75Figure 3.15 Performance comparison statistics over UV parameterizationproblems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Figure 3.16 Performance comparison statistics over 2D deformation prob-lems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78Figure 3.17 Performance comparison statistics over 3D deformation prob-lems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 4.1 A glass goblet broken into multiple shards example . . . . . . 84Figure 4.2 A graphical visualization of the rigid body fracture framework 89Figure 4.3 An illustration of boundary integral formulation . . . . . . . . 91Figure 4.4 Evaluation of different boundary integral discretization . . . . 94Figure 4.5 A 2D illustration of RBF based Fast Multipole Method . . . . 97Figure 4.6 An illustration of an RBF approximation to a 2D anisotropicexterior potential field . . . . . . . . . . . . . . . . . . . . . 98Figure 4.7 Runtime comparison between brute force method and FMM . 99Figure 4.8 Propagating fracture following a predefined path . . . . . . . 100Figure 4.9 Crack initiation local mesh operation . . . . . . . . . . . . . 102Figure 4.10 Crack extension local mesh operation . . . . . . . . . . . . . 103Figure 4.11 Crack splitting local mesh operation . . . . . . . . . . . . . . 103Figure 4.12 Crack merging local mesh operation . . . . . . . . . . . . . . 104Figure 4.13 Crack topology change local mesh operation . . . . . . . . . 105Figure 4.14 Three sculptures on a dinner table . . . . . . . . . . . . . . . 107Figure 4.15 Crack propagation inside a piece of transparent material . . . 107Figure 4.16 Rigid body brittle fracture simulation of plate and hollow sphere 109Figure A.1 Illustration of shape categories and topology operations . . . . 130Figure A.2 Multimesh data structure generation pipeline . . . . . . . . . 131Figure A.3 An illustrative example of multimesh structure . . . . . . . . 132Figure A.4 Comesh Optimization improves both the mapping and mesh-ing quality . . . . . . . . . . . . . . . . . . . . . . . . . . . 133Figure A.5 Boundary pairs construction pipeline illustration . . . . . . . 134xvFigure A.6 Illustration of MIPS energy’s transitivity property . . . . . . . 142Figure A.7 Example of blending dynamics into a flying bird animation . . 143Figure A.8 Evaluation of our animation inbetweening method . . . . . . 143Figure A.9 Evalution of our animation method by reproducing benchmarktasks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144Figure B.1 AKAP comparison to Newton-type methods . . . . . . . . . . 150Figure B.2 AKAP and SLIM comparison . . . . . . . . . . . . . . . . . 151Figure C.1 n particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160Figure C.2 quad tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161Figure C.3 well separated cell . . . . . . . . . . . . . . . . . . . . . . . 161Figure C.4 particle to multipole . . . . . . . . . . . . . . . . . . . . . . 162Figure C.5 multipole to multipole . . . . . . . . . . . . . . . . . . . . . 163Figure C.6 multipole to local . . . . . . . . . . . . . . . . . . . . . . . . 164Figure C.7 local to local . . . . . . . . . . . . . . . . . . . . . . . . . . 165Figure C.8 local to particle . . . . . . . . . . . . . . . . . . . . . . . . . 166Figure C.9 near field summation . . . . . . . . . . . . . . . . . . . . . . 166Figure C.10 particles with depth-6 quad tree . . . . . . . . . . . . . . . . 169Figure C.11 rank deficient interaction between well-separated clusters . . . 171Figure C.12 full rank pattern (colored in red) of Figure C.10 from tree level2 to 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172Figure C.13 low rank approximation of A represented in FMM . . . . . . 172Figure C.14 algebraic expansion of FMM . . . . . . . . . . . . . . . . . . 173xviAcknowledgmentsIt has been a great pleasure to spend almost six years in Vancouver studying andworking on some interesting and challenging math problems. I would like to thankDr. Robert Bridson for providing me enough freedom to explore whatever I aminterested in as well as your patience and support when I was struggling. I also wantto thank Dr. Chen Greif for agreeing to be my co-supervisor to keep me in track.I am really appreciated the discussion and brainstorm we had in your office onnumerical problems. I will also thank Dr. Essex Edwards. I cannot imagine what itwould be like if you were not in the lab when I started the journey. Thanks a lot forbeing patient when I bothered you with so many questions everyday. Thanks alsogo to Dr. Danny Kaufman and Dr. Jovan Popovic. Thank you guys for providingme the internship opportunity at Adobe, where I not only had two published papersbut also met many talented people and friends. I want to thank Dr. Xinxin Zhangfor being a great lab mate as well as room mate. I enjoyed learning fluid mechanicsas well as visual effects from you. Moreover, I would like to thank my friends Dr.Lei Xiao, Dr. Yifan Peng, Dr. Jianhui Chen, Dr. Lili Meng, Dr. Qingkun Su, Dr.Xin Tao, Dr. Chengzhou Tang, Chi Zhang, Dr. Jason Peng, Ye Fan, Dr. QiangZhou, Dr. Edwin Chen, Dr. I-Chao Shen, Dr. Minchen Li, Dr. Shuochen Su,Weicong Liao, Dr. Rui Ge, Dr. Chunsheng Zhu, Yajian Tong, Dr. Shan Jiang, Dr.Guang Gao, Dr. Todd Keeler, Dr. Chongyang Ma, Dr. Libin Liu, Dr. Will Chang,Dr. Wei Sun, Dr. Hong Hu. Thank you guys for sharing your research work withme. It’s always fun to learn new stuff from people who work in different areas.Finally, I want to thank my parents. Thank you very much for all your support sofar. 感谢爸妈！xviiChapter 1IntroductionThe mapping problem is, given geometric inputs and a possibly under-specifiedoutput (such as just its boundary, or even just a restriction to the plane), to find a“good” function or mapping that takes the inputs to the outputs. It plays a fun-damental role in the computer graphics field. After introducing specific measure-ments of what “good” might mean, researchers generally look for an optimal so-lution, potentially also subject to some problem constraints. For example, the UVparameterization problem is to find an optimal mapping between a surface mani-fold embedded in the 3D space and a 2D plane, with notions of optimality that mayinclude area preservation, angle preservation, or distance preservation. A locallyoptimal solution can be obtained by minimizing the mesh distortion after flatten-ing the surface mesh into the plane as illustrated in Figure 1.1. Other interestingand popular applications include planar animation and keyframe interpolation (thesubject of the first section of this thesis), image warping, and 3D modeling andelastic simulation (the subject of the last section of this thesis, with the additionalcomplexity of topology change in the form of fracture). Such problems usually canbe formulated as a (potentially nonlinear) optimization and may require extensivecomputation to achieve acceptable results. The situation can be even more chal-lenging when topology change must occur during the optimization. For example,finding an optimal mapping between topologically in-equivalent shapes, as shownin Figure 1.2, is non-trivial, as a singularity problem needs to be resolved before acontinuous map can be established.1Figure 1.1: A UV parameterization example: a curved surface mesh forms ahand geometry in 3D and is cut into a disk topology. A local optimiza-tion finds a 2D planar embedding with minimal isometric distortion,shown on the right.We start our discussion on optimal mapping problem through a 2D animationsystem introduced in Chapter 2. We propose a variational shape interpolation basedapproach to reduce the amount of manual work in the existing animation pipeline.To construct the interpolation path that represents the minimum amount of workrequired to transform from one keyframe shape to the other, an optimal mappingneeds to be established between keyframes first. To account for topology differ-ences between keyframes, and hence topology change during interpolation, wepropose a novel multimesh data structure that resolves the topology inconsistencyfrom the given keyframe shapes. Then a Comesh Optimization framework is intro-?Figure 1.2: There is no bi-continuous map between shapes with differenttopology. To establish a mapping between such topologically inequiva-lent shapes, singularity problems need to be resolved as highlighted inthe dotted region.2duced to build an optimal map in terms of both intra-mesh and inter-map quality.This framework is composed of two nested loops: the outer loop walks through theinput shapes one by one and optimizes the geometry and topology of each shapealternately; in the inner loop, we adopt a simple gradient descent method to op-timize the mesh geometry according to an appropriate energy function measuringdistortion.The two loop optimization framework provides satisfying results in the end,but suffers from slow convergence due to the adoption of gradient descent methodin the inner loop. Therefore we turn to nonlinear optimization techniques in Chap-ter 3, proposing a new solver with state-of-the-art performance. Our solver notonly shows thoroughly superior or at least competitive performance across a largerange of mesh distortion optimization problems, but also resolves some lingeringrobustness problems in prior solvers related to termination criteria, with a reliableand largely mesh- and energy-independent convergence test.One limitation of our solver is that, in 3D, its performance scales with thevolumetric mesh resolution, even for applications where only the geometry of thesurface matters ultimately, such as in the physics-based animation of elastic solids.If we include fracture in the simulation, volumetric meshes further present seriousperformance difficulties as nontrivial conforming tetrahedral remeshing is requiredaround dynamic cracks. To account for these two issues, we propose a boundaryintegral technique in Chapter 4, formulating the optimization problem underly-ing quasistatic elasticity over the boundary mesh (rather than the interior volume),eliminating the need for a volumetric mesh and with it the related performancescaling trouble and volumetric meshing problem. We also introduce a kernel-independent Fast Multipole Method to solve the new formulation efficiently. Bycombining this approach with an explicit surface-tracking framework, topologychange of the simulation domain can be handled in a more elegant manner. Weshow such approach is especially useful when applied to brittle fracture modeling.We describe the major technical details in the following three chapters. Somederivation and proof details are provided in the appendix.3Chapter 2Planar Interpolation withExtreme Deformation, TopologyChange and Dynamics2.1 IntroductionHumans draw what’s on their mind: they do it to tell stories; they do it to sharewith others; and they do it to reflect on their thoughts. They use cameras whenthey record, but they draw when they invent. The software industry has long rec-ognized this need with the development of comprehensive products for sketching,illustration, and presentation. Brushes, layers, undo, copy, paste, and other toolscan make it easier to draw. However, when it comes to drawing animation, eventhe most sophisticated products provide little assistance for changing shape. Inthis work, we develop a new method to enable artist-driven inbetweening of drawnshapes as an essential steps towards drawing animation without having to drawevery frame.Animations change shape. Drawing software requires redrawing each frame.Animation software does not; but it requires rigging to articulate shape, and keyfram-ing to change that articulation over time. A drawer wants to draw - not rig orkeyframe - but not every frame. This gap inhibits animation. Almost every child4Figure 2.1: We interpolate inbetweens (unboxed images) from sparse keydrawing shapes (boxed images) across arbitrary topology changes andextreme deformations. Artists define desired correspondences to ex-plore interpolation paths in-between key shapes and add dynamic effectswhen desired (bottom).draws pictures, but almost none animate.Making animation more accessible today does not require eliminating draw-ing. It requires assistance with the creation of inbetween drawings. We focus onenabling intuitive, rig-free control of inbetweening from a set of drawn shapes. Oursystem delivers an initial animation from just a few sparse key drawings and a setof user-specified correspondences, see Figure 2.1 and Section 2.3, and then can berefined with swapped-in and added drawings, e.g., extreme or breakdown poses.This emulates the well established workflow from classical (hand-drawn) anima-tion, but with our software assisting with the interpolation of inbetween drawingshapes.Inbetweening with drawings should allow artists to communicate as freely asthey draw. Drawings can exaggerate and invent without constraint. Thus com-putational inbetweening with drawings should likewise be unrestricted by riggingand articulation curve constraints. Current computational tools, however, cannotdeliver this today because geometric methods do not yet support this implied fullfreedom of artistic expression. Even the most simple outline drawing shapes revealthree difficult-to-inbetween transformations: 1. large and exaggerated non-uniform5stretching and compression; 2. topological change and 3. physics-based dynamics.Contributions To enable artistic freedom of expression, we remove these priorrestrictions on the range of possible changes between shapes supported by compu-tational interpolation methods. We present a set of contributions that support mesh-based, variational interpolation for interactively creating directed inbetweens fromarbitrary and unrestricted sets of 2D drawing shapes without rigging. Our goalthen is to deliver pleasing interpolation that is smooth, stable, and consistent withextreme deformation, unrestricted topology change, and physics-based dynamic ef-fects. More specifically, if the shape interpolation path is parameterized in someway (detailed in later chapter), stability and consistency are defined in the senseof continuously depending upon time and the input geometry (no abrupt changesfrom small perturbations), where smoothness is refered as being differentiable withrespect to the parameterization.Solving this problem led us to two mutually dependent subproblems: 1. aninterpolation method for pleasing transitions requires higher quality compatiblemeshes than previously available and new energies to support these transitions;and 2. compatible meshing must be aware of the interpolation method applied.Our primary contributions address these two subproblems:• a Comesh Optimization algorithm that improves both the mesh and mappingquality of compatible meshes so that we can support pleasing interpolationacross unrestricted shape changes; and• a Multimesh Interpolation method that extends variational interpolation [111,137] by introducing a consistent, annotated multimesh structure and a newshape-space energy that supports topology change, can prevent shape overlapwhen and where desired, and add dynamic effects as directed. The result-ing interpolation is efficient and enables interactive animation creation andediting. We demonstrate this with a range of challenging animation tasksincluding production examples.62.2 Related Work2.2.1 Planar AnimationPlanar animation is classical topic in computer graphics. With decades of explo-ration in the area, a wide range of methods have been developed for 2D animationusing strokes [67, 142], triangle meshes [2, 10, 11] and vector graphics [36], toname just a few domains. To enhance planar animation, researchers have alsoconsidered occlusion in 2D. Yang et al. [152] proposed a multi-element 2D shaperepresentation to resolve disjoined and overlapped parts. Sy´kora et al. [133] intro-duce cracks in the embedding cage by allowing discontinuities in the depth field.Embedding cages can then be separated in two in order to improve image regis-tration accuracy. Nevertheless, these approaches cannot model topology changessuch as adhesion effects during interpolated animation. Whited et al.[142] resolvechallenges in layering and support occlusion ambiguities with the help of artists.Research in planar animation has also explored velocity field advection [74] anddata driven strategies [6], inference of motion cycles from motion snapshots ofindividuals captured in a stills [151], and beautification of existing drawings viaglobal similarity analysis [149]. Here we introduce a novel interpolation methodto enable planar animations from sparse key shapes with extreme deformation,topology change and dynamics—to our knowledge this is the first time this hasbeen possible.2.2.2 Shape InterpolationShape interpolation has been widely explored in geometric modeling, computer vi-sion, imaging and computer animation [146]. As in Euclidean space, the action ofinterpolation on shapes is induced by definition of a metric such as Minkowski [62],Wasserstein [122], and Gromov-Hausdorff [45] distances, as well as local distor-tion measures [48, 146]. For design purposes, interpolation should satisfy artisticintent. Yet, as a geometric task, interpolation between arbitrary planar shapes isambiguous and underspecified. Consider, for example, the classic exercise of ani-mating a square to an arch. Depending on taste and storyline, we may either wishto stretch the square out in anticipation, bridging to the arch (Figure 2.2 stretch);7StretchSplitFigure 2.2: Animating a square to an arch: depending on taste and storyline,an animator may wish to (1) stretch the cube out in anticipation, bridg-ing to the arch; or (2) may instead want to first split the cube and thenspread to the arch. We provide flexible exploration of the range ratherthan being restricted to one animation pre-established by a metric as ine.g., optimal transport which captures the split but misses the stretch.or we may instead wish to first split the square from the bottom and then spread tothe arch (Figure 2.2 split). Our method follows the local distortion framework, aswe seek a method that provides control via correspondences over interpolation andtopology changes so that users can control interpolation path.Previous work includes comprehensive review of local distortion approaches[28, 137], but we prefer a categorization into 1. geometric and 2. physical ap-proaches. Geometric approaches rely on purely geometric quantities such as thepullback metric tensor (i.e., rotation-invariant right Cauchy-Green strain); Greenstrain or edge length; dihedral angles; or a geometric measures of shape distortion[28, 65, 72, 81, 144]. Our approach follows variational shape interpolation, basedon physical quantities [26, 48, 49, 111, 137, 145, 146]. Unlike previous approachesin this family, our final method does not require a full correspondence map fromusers and supports interpolation with topology changes.82.2.3 Compatible MeshingAs in traditional mesh-based interpolation, our method requires a compatible mesh-ing of all example shapes. Shapes with moderate deformation require no more thanexisting shape modeling techniques [8]. However, shapes with large and exagger-ated deformations and/or topological changes need higher-quality meshes adaptedfor these changes. We introduce new methods here to efficiently generate high-quality compatible meshes for collections of 2D shapes, given only a sparse set ofcorrespondences.Locally injective mapping has been extensively studied [11, 75, 104, 113, 132,139]. Our problem setting is most similar to those investigated by Weber andZorin [139]. Here, however, we focus only on planar shapes arising from draw-ings. To do so we must solve a problem currently unaddressed by previous meth-ods: guaranteeing local injectivity is necessary but not sufficient to obtain stable,consistent and smooth interpolations over extreme shape changes. Locally injectiveinitializations require additional optimization for both mesh quality and mappingquality.Prior methods on forming compatible meshes [11, 71, 106, 107, 112, 128–132]similarly observe the need for optimizing meshes after initial connectivity is estab-lished. However, these prior methods optimize only the quality of each mesh anddo not fix mapping distortions between them. Unfortunately, while optimizationof just mesh quality can occasionally be sufficient for small deformation interpola-tions, it is generally insufficient, see e.g., Figure 2.8, especially for large deforma-tion interpolations where mapping distortion is a critical issue. For this purpose wemotivate and introduce Comesh Optimization a new mesh optimization process forhigh-quality compatible meshing that improves both mesh and mapping quality.2.3 Planar InterpolationAs shown in Figure 2.3, we begin with a set of drawings, i ∈ [1,m], in the planeas input. We then provide an interactive tool to extract drawing shapes and tospecify desired correspondences for cuts, openings and boundaries. We next applyour Comesh Optimization algorithm to compute compatible, high-quality meshes,each with n vertices Xi ∈ R2n, that deliver pleasing interpolations (smooth, stable,9Input Drawings Multimesh Constructiona bcdea a abb bcc cdddeeedCut OperationCorrespondenceecComesh Optimization Animating InbetweensFigure 2.3: Our framework takes drawings as a set of images. It constructsour multimesh structure by pairing extracted shapes with annotated cutsand correspondences to build initial embedding shapes for each withcompatible meshing. Our Comesh Optimization then improves bothmesh and mapping quality of the shape to obtain high-quality, stableinterpolations.and consistent) when minimizing our interpolation energy. With these generatedmeshes in hand, our framework provides user controls to interactively explore in-between synthesis by interpolating with varied timings, energy parameters, inter-polation weights and artwork sequencings.2.3.1 Inbetweening via Variational InterpolationWe compute inbetweens by interpolating variationally [26, 48, 111, 145] betweenour m example artwork shapes. We start by defining a potential W (X ,x) to measurethe energy of deforming shape x away from a reference shape X . Variational in-terpolation between just two shapes, X1 and X2, follows the approximate geodesicpath [26, 48] between them asx¯(w) = argminx(1−w) W (X1,x)+w W (X2,x), (2.1)with w ∈ [0,1].Variational interpolation over many shapes is extended by shape averaging [137,145]. Constructing a weighted sum of deformation energies,WI(x,w) =m∑i=1wi W (Xi,x), (2.2)10we interpolate with the minimizationx¯(w) = argminxWI(x,w), (2.3)varying convex weights w = (w1, ..,wm)T ∈ Rm+, ∑mi=1 wi = 1. This multishapeinterpolant is then effectively an equilibrating deformation that balances betweenweighted contributions of rest shapes. As weights w vary, we move between theinfluence of multiple shapes to create blended inbetweens.Inbetweeing For a set of example artworks A = {X1, ..,Xm}, a single inbetweenis given for each choice of weights w ∈ Rm+. To build animations, we then candesign sequences of inbetweens with animation curves w(s) ∈ Rm+, s ∈ R appliedin x¯(w(s)). Keeping just two nonzero entries in w at a time retrieves a traditionalkeyframe inbetweening, while blending across many nonzero weights at a timeestablishes complex inbetweens from multiple artworks simultaneously. See Sec-tion 2.6 for details.Extending Variational Interpolation Variational interpolation promises a power-ful approach for freeform inbetweening with many shapes. However, methodsbased on variational interpolation have so far been generally restricted to shapespaces with a limited range of expression. Shapes used in these interpolationscould neither contradict one another, nor strongly counteract the underlying de-formation energy [81]; e.g., by extreme deformation, change in topology, or evenscaling, see Figure 2.4 for examples.We address these limitations to enable free-form interpolation between shapesin three steps. First, we introduce here a multimesh structure that ensures a com-patible meshing across all shapes and enables annotation of desired correspon-dences for opening (respectively closing) and cutting (respectively merging) be-tween shapes. Second, we augment the standard elastic measure of deformation inthe shape average with additional energies that measure change on boundary: sep-aration/merging and overlap. When minimized in the shape average, these ener-gies now allow interpolation over topology change and, when desired, additionallypreserve non-overlap in artworks. Third, we propose a compatible meshing algo-11Unoptimized Interpolation Comesh Optimized InterpolationSquare - SquarePentagon - Slim StarPentagon - StarFigure 2.4: Interpolated inbetweens generated by LIM parameterization (leftcolumn) and Comesh Optimization (right column). As we explore shapevariation, e.g., thickness of the star geometry in the two pentagon-to-starexamples above, interpolation behavior should vary consistently. With-out Comesh Optimization (left) interpolating with an initial LIM pa-rameterization generates unsatisfying animation artifacts, non-smoothinterpolation sequences, and inconsistent results for inbetweening. Af-ter Comesh Optimization (right) both mesh and mapping quality are im-proved - we obtain low-distortion, smooth and consistent interpolationover varying input shapes.rithm, Comesh Optimization, that critically improves the quality of our multimeshto obtain high-quality, artifact-free interpolations over extreme deformations andtopological changes.2.3.2 MultimeshWe build our multimesh structure with a single shared mesh topologyO = {V,E,T}of vertices, V , edges, E, and faces, T , and m differing vertex geometries, Xi ∈ R2n,cohesive edge sets, Ci ⊂ E×E, and boundary vertex sets, Bi ⊆ V for all exampleshapes i ∈ [1,m]. Thus, for example, vertex vk ∈ V will have different positionsXki = Xi(vk) and Xkj = X j(vk) ∈ R2 in shapes i and j respectively.To enable topology change we introduce cohesive edge sets, Ci, per exam-ple shape. These are sets of cohesive edges, C ji : edge-edge pairs on the interior(respectively boundary) of each example shape that indicate edge-pair correspon-dences where a shape will open (respectively close) or cut (respectively merge) forinterpolation between topologically inequivalent shapes, see Figure 2.5.12~(a) cut operation~(b) open operationFigure 2.5: Non-self-intersecting planar shapes are converted to the sameequivalence class, with continuous deformations possible betweenthem, via two topological operations: (a) cutting (respectively merging)and (b) opening (respectively closing). Cut and opening locations areannotated by cohesive zones defined in shapes’ interior region. Upondiscretization, each cohesive zone is a set of cohesive edges.Currently, in our framework, cohesive edge sets are created by artists draw-ing their desired correspondences. As we will see this enables the exploration ofan expressive range of interpolation. Looking ahead we have also designed ourmultimesh structure and interpolation so that it will similarly function if corre-spondences are alternately discovered by other interactive or automatic workflows.Once designated, cohesive edges C ji = {{va,vb},{vc,vd}}match generally non-adjacent edges (Xai ,Xbi ) and (Xci ,Xdi ) for correspondence in shape i; see Figure 2.6.These cohesive-zone edge-pairs are then pulled towards and separated away fromone another in interpolation via a cohesive energy we construct in Section 2.5. Foreach shape i, we then likewise gather all vertices on its boundary into a boundaryvertex set Bi ⊆V ; we use these vertices to prevent artwork overlap, when desired,during interpolation.We then have multiple compatible boundaries of non-simply connected planarshapes (multiply connected and disjoint). Compatible piecewise linear maps, i.e.triangle meshes, could then be intuitively established in three ways: 1. triangulateone shape and map it to the others [75, 113]; 2. compatibly triangulate multipleshapes simultaneously [11, 132]; or 3. triangulate each shape separately and thenreparameterize them compatibly [139].Here we adopt the first approach and use Triangle [114] to generate a con-forming Delaunay triangulation of one shape first. As this meshing step may in-sert Steiner points along the boundaries, we then upsample boundaries of the other13Cohesive EdgeFigure 2.6: A cohesive edge, C ji , is a pair of collocated edges which can sep-arate apart under deformation. The induced displacement field in Xi willbe discontinuous over C ji .shapes and keep any splits in cohesive edges consistent. Finally, we adopt Schuelleret al.’s LIM algorithm [113] to establish locally injective mappings of the trianglemesh to all other shapes. We note that this approach does not guarantee an ini-tializer for arbitrary cases; we could, for instance, follow the third approach asour workflow can start with results produced by any of the above three strategies.However, in our experiments and for our examples, we have so far found strategy 1with LIM to be most effective to bootstrap our multimesh construction. For furtherdetails on the creation of our multimesh structure with correspondences we referthe interested reader to Appendix A.In Section 2.5 we will construct and detail our deformation energy over meshes,cohesive energy over cohesive edges, and non-overlap energy over boundary ver-tices. However, the resulting interpolation we will obtain using the default multi-mesh will lead to unsatisfying animation artifacts during interpolation (see e.g.,Figure 2.4) due to, as we will soon see, both poor intra-mesh element qualityand poor inter-mesh mapping quality. In Section 2.4, we thus first introduce ourComesh Optimization algorithm to improve both the mesh and mapping qual-ity of our multimesh structure to resolve these interpolation artifact problems forfreeform shape collections. Here and in the following sections we adopt the nota-tional conventions summarized in Table 2.1.14Table 2.1: NotationNotation DefinitionV shared mesh vertex setE shared mesh edge setT shared mesh face setn = |V | number of verticesO = {V,E,T} shared mesh topologyXi ∈ R2n ith example-mesh’s vertex positionsm number of example shapesA = {X1, ..,Xm} set of example-shape meshesx ∈ R2n deformed-mesh vertex positionsw ∈ Rm+ interpolation weightss ∈ R timing parameterw(·) : R→ Rm+ animation curvex¯(·) : Rm+→ R2n variational interpolationak(·) : R2n→ R signed area of triangle tk ∈ O in input meshCi ⊂ E×E cohesive edge set for ith meshBi ⊆V boundary vertex set for ith meshX ji = Xi(v j) ∈ R2 vertex v j ∈V in mesh Xix j = x(v j) ∈ R2 vertex v j ∈V in deformed mesh x2.4 Comesh OptimizationArtifact-free interpolation over changing shapes requires high-quality mappingsbetween well-formed, compatible meshes. Mesh quality, per shape, is essentialfor the stability of minimizations performed on discretized energies for variationalinterpolation. Mapping quality, between shapes, ensures local similarity for high-quality visual interpolations. In this section we propose a Comesh Optimizationmethod that jointly improves both mesh and mapping quality of initial, inversion-free multimeshes. Our resulting optimized multimeshes then give smooth, stableand consistent interpolations that are visually pleasing over extreme changes inshape.Goal Given multiple, compatibly parameterized boundaries of non-overlappingplanar shapes we seek meshes with 1. individually well formed elements and 2.bi-directionally optimal piecewise-linear maps between all shape pairs.15Figure 2.7: A low resolution example highlights how initial, locally injective,compatible meshes of example shapes (pentagon and star on far left andright) will produce poor interpolation results (top); note sliver trianglesand large distortion maps between triangles. Comesh Optimization ofthis initial mesh then generates pleasing interpolations (bottom) afterimproving triangle quality and mapping distortion of the meshes.Input We bootstrap this process, as described above in Section 2.3, with an ini-tial, locally injective parameterization of shapes annotated with correspondences.We apply Schu¨ller et al’s [113] LIM algorithm to obtain a single mesh topology,O ,for m differing shapes with varying vertex positions, Xi, i ∈ [1,m]. This commonconnectivity delivers dense mappings from the specified correspondences definedin Section 2.3.2.Other options for bootstrapping to an initial connectivity are also possible [11,132, 139]. However, while initial meshes found by these alternate approaches arealso often locally injective, they too have no guarantee of quality. In practice out-put meshes from all such methods, including LIM, require additional optimizationto avoid interpolation artifacts; see e.g., Figures 2.4, 2.7 and 2.8.We have so far found LIM to most consistently deliver good initializations forour approach. Neverthelss, on its own, LIM may produce unacceptable artifactsupon interpolation if the energy is inappropriate, or the chosen direction of themapping causes troubles, see e.g., Figures 2.4 and 2.7. This motivates our ComeshOptimization algorithm for improving the quality of both mesh element shapes and16Before AfterFigure 2.8: Optimizing only the mesh quality of example shapes and ignor-ing map distortion as in [11, 132] improves triangle shape but leavesdistortion (colormap) largely unimproved (top right). In turn, resultingvariational interpolations in between the shapes (bottom row) then suf-fer from symmetry breaking, inconsistency and instabilities. Comparewith our Comesh Optimized results in Figure 2.4.inter-mesh mapping.2.4.1 Mesh and Mapping EnergiesWe measure the mesh quality of each shape i’s triangulation with Alliez et al.’s [3]Optimal Delaunay Triangulation (ODT) energy,D(Xi,O) =14 ∑v j∈V‖X ji ‖2Γ j(Xi)−qi, (2.4)and measure the map distortion from shape i to shape j with the Most-IsometricParameterizations (MIPS) energy [54],E (Xi,X j,O) = ∑tk∈Tak(Xi)(σ1kσ2k+σ2kσ1k). (2.5)17For the variational Delaunay ODT energy above, Γ j(Xi) is the one-ring area ofvertex X ji while qi is a constant term for fixed-shape boundaries1. For the MIPSenergy, ak(Xi) is the area of triangle tk, in shape Xi, while σ1k ,σ2k are the first andsecond singular values of the 2×2 deformation gradient for triangle tk ∈ T , whentreating shape X j as a deformation away from shape Xi.2.4.2 Optimizing Solely Mesh or Mapping is InsufficientComesh Optimization is built from our observations that 1. only optimizing trian-gle quality per mesh, as done previously, will not prevent high-distortion mappingsbetween meshes; and 2. only minimizing mapping distortion between meshes,irrespective of choice of distortion measure, will not guarantee good-quality trian-gulations for stable interpolation.Minimizing solely distortion, irrespective of distortion measure, results in am-biguity: many mesh solutions are possible and they rarely provide well-shapedtriangles; see below. We choose to jointly apply the ODT energy and a cyclic sum(detailed below) of MIPS energies to optimize triangle quality and map distortionrespectively. For measuring distortion in Comesh Optimization we initially investi-gated other energies including the isometric energy from Smith and Schaefer [119]and ARAP [124]. However, they both resist scaling and so fight the exaggeratedscaling behaviors we want to allow. In the end we choose conformal MIPS as itlimits angular distortion and allows large scale changes while preventing inversion.For very similar shapes, e.g., between squares of similar orientation and size,if we optimize just mesh quality alone, we will sometimes also obtain a close-to-optimal, low-distortion map. However, for a general collection of different shapes,optimizing solely mesh quality will produce a set of meshes that are individuallywell formed yet work poorly for interpolation due to large map distortions betweenthem. See for example Figure 2.13 (ODT), where we optimize solely mesh qual-ity with the ODT energy and so still obtain high distortion; correspondingly, inFigure 2.8 we see that optimizing just mesh quality introduces unsightly artifactsduring interpolation.On the other hand, if we optimize only map quality, e.g., here with the MIPS1In our setting qi thus has no effect over gradients with respect to interior vertices, but its inclusiondoes cancel out gradient components due to boundary vertices, see Alliez et al. [3].18Similar Input ShapesConformal Map 1Conformal Map 2Figure 2.9: Ambiguity in optimal, low-distortion, piecewise-linear maps:here both maps minimize the same distortion measure (in this caseMIPS) yet Map 1 has better mesh quality.energy, our problem is ill defined. Consider Figure 2.9 where both maps are con-formal yet they present different meshes. There exists ambiguity in optimally low-distortion piecewise-linear maps between shapes and so, unless we are exceed-ingly lucky, optimizing for map quality alone will not give us well-formed meshesfor pleasing interpolation. For example, in Figure 2.13 (MIPS) optimizing withthe MIPS energy we obtain poor-quality meshes even though we generate a low-distortion map.These issues only worsen as we consider general, multishape interpolationsbetween m > 2 shapes and so we design our Comesh Optimization method tooptimize both mesh quality with ODT and map quality with MIPS. However, tooptimize map quality we must somehow efficiently minimize mapping distortionbi-directionally between all shape pairs from our m example shapes. We addressthis challenge next.2.4.3 Multishape MIPS and SymmetryAs we have many shapes to be compatibly meshed, our goal is distinct from manytraditional mesh optimizations task. Each example shape is effectively treated as arest shape for measuring distortion energies in variational interpolation and so re-quires well-formed mesh elements. At the same time, for measuring mapping qual-ity, each mesh is also a deformed shape with respect to all other example shapes inour collection.In our setting interpolation is thus not a one-way trip from one shape to an-other. It can take any direction in shape space with any combination of shapesparticipating and so our mappings should likewise minimize distortion for all pos-19sible paths between shapes. Consider again the simple square-to-hexagon examplein Figure 2.13 where square and hexagon respectively have vertex positions Xsqand Xhe. Following standard distortion minimization we first try optimizing overour target hexagon’s vertices, minXhe E (Xsq,Xhe,O). This, however, gives us the“Asym” solution in Figure 2.13 where we see that distortion after this optimizationremains poor. Why does this happen? The full degrees of freedom in this problemhave been ignored: we can and should optimize distortion in all directions by treat-ing each shape as both a rest and deformed shape when we optimize total mappingquality.Cyclic Summation Hence we seek symmetric, optimally low distortion maps forboth Xi→X j and X j→Xi ∀i, j ∈ [1,m]. For square-to-hexagon this simply amountsto minimizing the sum of E (Xsq,Xhe,O) and E (Xhe,Xsq,O) giving us the low dis-tortion solution in the “MIPS” row of Figure 2.13.More generally we seek to minimize all possible distortions between examplepairs. While this is straightforward to do for two example shapes as above, opti-mizing this explicitly as an objective over all pairwise distortion measures betweenm example shapes is prohibitively expensive.We side-step this potential computational obstruction by observing that MIPSis transitive: if the conformal distortion between Xi and X j and between X j and Xkare both small (bounded above), then the conformal distortion between Xi and Xkis also small (bounded above). See Appendix A for proof and details. This allowsus to construct our symmetric, all-pairs distortion energy with just a cyclic sum ofthe MIPS energy over our m example shapes asm∑cycE (Xi,X j,O) =E (X1,X2,O)+E (X2,X3,O)+· · ·+E (Xm−1,Xm,O)+E (Xm,X1,O).(2.6)In practice we observe stable traversals of the shape space when optimizing withour cyclic MIPS energy, ill-posed results without it, and also confirm in experimentthat we obtain the same interpolation behavior even as we vary the order of exampleshapes in the cyclic summation.20Input Output50 100 15050010001500Example Cell Splitting50 100 15050010001500200025003000Angle Distribution26+4Figure 2.10: Comesh optimization improves both map and mesh quality ofour multimeshes subject to position and topology constraints. Herewe visualize the change in map quality (colormap) and mesh quality(angle distribution) for our cell splitting example. Designated cohesiveedges (red) are preserved by both vertex improvement and topologyupdate steps.2.4.4 Comeshing ConstraintsWe have so far focused on constructing a suitable objective for optimizing our mul-timesh. To maintain prescribed correspondences, i.e., points and cohesive zones,and to maintain consistent mesh topology we must also impose constraints on ouroptimization. We begin by constraining feature points on boundaries and corre-spondences that can be specified as positional constraints on vertex positions pershape, P(·) = 0. To preserve multimesh structure we implicitly constrain meshtopology O to remain consistent across all meshes. We then enforce local injec-tivity, ak(·)> 0, explicitly so that topology optimization steps (e.g., edge flips andedge collapses, see below) do not cause our inter-mesh mapping energies, Equa-tion 2.6, to go towards infinity as area becomes small.Finally, to maintain pre-defined cohesive zones, we constrain edge topologyand vertex positions to stay equivalent on both sides of all cohesive zones. Recallthat each cohesive edge C ji corresponds to two edges in the edge set E, with respec-tive end vertices, {Xai ,Xbi } and {Xci ,Xdi } in shape i. We constrain all such pairs tobe collocated, during all topological steps with the constraintZ (·, ·) = 0, requiringXai = Xci , Xbi = Xdi . (2.7)See Figures 2.10, 2.11 and 2.12.21nj-th cohesive edgeFigure 2.11: We require the two pairs of end vertices on a cohesive edge tobe collocated.2.4.5 Comesh OptimizationIn summary our Comesh Optimization seeks a piecewise-linear map that jointlyminimizes our cyclic MIPS and area-scaled ODT energies to find a compatible De-launay triangulation satisfying position, cohesive zone, and local injectivity con-straints,minX1,X2,··· ,Xm,Om∑i=1D(Xi,O)∑tk∈T ak(Xi)+m∑cycE (Xi,X j,O)s.t. P(Xi) = 0, ∀i ∈ [1,m],ak(Xi)> 0, ∀tk ∈ T and i ∈ [1,m],Z (Ci,Xi) = 0, ∀i ∈ [1,m].(2.8)Here note that we normalize each shape i’s ODT energy by its area ∑tk∈T ak(Xi).This scales the energies D and E in our objective to common units of squaredlength.A note on relative weighting Comesh Optimization minimizes a multi-objectivefunction formed by the sum of the cyclic MIPS energies for mapping and the area-scaled ODT energies for meshing. Note that we choose to scale both terms equallyin our objective above in Equation 2.8. We base this choice on our prior observationthat the space of conformal piece-wise linear mappings is generally ambiguouswith respect to mesh quality. This, in turn, suggests that we can find a wide range22Inverted TriangleFlip BoundaryFlip Cohesive EdgeInverted TriangleDegenerated BoundaryModi!ed BoundaryEdge FlipEdge SplitEdge CollapseConsistent DiscretizationInconsistent Discretizationedgecohesive edgeedge to be modi!edvertexvertex to be added reference meshTopology Update Cohesive Edge Discretization Invalid Edge Flip Invalid Edge CollapseFigure 2.12: We apply three local operations to update mesh topology: edgeflips, edge splits and edge collapses. These operations are constrainedin order to satisfy our Comesh Optimization constraints. Each opera-tion is applied to the multimesh topology shared by all example shapemeshes. It propagates to all meshes (or none) to maintain consistenttopology across all shapes.of conformal mappings with differing mesh qualities. To better understand howthe relative weighting of these two terms, mapping and meshing, can actually varyin practice in our results, we experimented with changing their relative weightsin our optimizations. For all combinations of finite, nonzero weights we observequalitatively and quantitatively, see Appendix A, similar results from our ComeshOptimization and so keep them equally weighted in our objective.2.4.6 Solving Comesh OptimizationTo solve Comesh Optimization we construct an algorithm that minimizes our con-strained problem Equation 2.8 with block, cyclic descent over example shapes[1,m]. Our algorithm is detailed in Algorithm 1. Our outer loop iterates block-wise, by shape. For each outer iterate on shape i ∈ [1,m] we solve for an update ofshape i’s vertex positions Xi and shared multimesh topology O , while holding allother X j, j 6= i fixed, subject to our compatibility and feasibility constraints. Then,within each inner iteration on shape i, we optimize mesh Xi and topology O apply-ing one step of vertex optimization with gradient descent, and three steps of meshoptimization and refinement applying edge flips, edge splits and edge collapses.All four operations are constructed to preserve local mesh constraints and maintain23Algorithm 1 Comesh Optimization1: Input: multimesh data structure O,{Xi}, {Ci}, {Bi}, i ∈ {1, · · · ,m}.2: output: O,{Xi}, {Ci}, {Bi}.3: Initialize: max outer itr = 500, max inner itr = 100,4: outer tol = 10−6, inner tol = 10−7,5: [·] : Z→ Z/mZ.6: for loop outer = 0 to max outer itr7: total res = 08: for i = 1 to m9: inner res = 010: for loop inner = 0 to max inner itr11: edge flip(Xi, O , Ci, Bi)12: edge split(Xi, O , Ci, Bi)13: sub res= gradient descent step(X[i−1], Xi, X[i+1])14: edge collapse(Xi, O , Ci, Bi)15: if inner res< sub res16: inner res = sub res17: end if18: if sub res< inner tol break19: end for20: total res = total res + inner res21: end for22: if total res< outer tol break23: end forcompatible meshing across all m shapes.Each vertex optimization step improves current mesh coordinates by applying asingle step of projected gradient descent on vertices Xi. The gradient of our objec-tive in Equation 2.8 with respect to shape, Xi, is ∇Xi [D(Xi,O)∑tk∈T ak(Xi)+E (Xi,Xi+1,O)+E (Xi−1,Xi,O)]. For a detailed derivation of this gradient see our supplement. Eachgradient descent step applies bisection correction, per vertex, to find a conserva-tive improvement along the gradient that both respects our constraints on vertexposition and does not invert triangles in the one-ring neighborhood.Mesh optimization and refinement applies edge flips, edge splits and edge col-lapses. Edge flips improve mesh connectivity. For each candidate edge, we accepta flip when it both decreases the Delaunay objective over all meshes and respects24Input:CO:6+4MIPS:ODT:Example Square-HexagonAsym:50 100 15020406050 100 1502040608010012050 100 1502040608050 100 15020406080100120Angle DistributionFigure 2.13: Here we visualize the change in map quality (colormap) andmesh quality (angle distribution) for a range of multimesh optimiza-tion alternatives. Minimizing only MIPS accounts for only distortionand so still obtains poor quality elements (MIPS) while a solution min-imizing ODT considers solely mesh quality and so obtains poor map-ping quality (ODT). Likewise, minimizing distortion asymmetrically,in this case from just square to hexagon, will still produce poor so-lutions as in (Asym); compare to the (MIPS) solution that minimizesdistortion symmetrically. Our Comesh Optimization seeks a symmet-ric solution optimizing both per-shape mesh quality and all shape-pairmappings and so obtains both low distortion and well-formed triangleelements (CO).25our constraints in Equation 2.8. We alternate between evaluating feasibility andevaluating energy decrease. For our feasibility evaluation we check 1. if the flipcandidate is a boundary or cohesive edge; and 2. if flipping leads to an inversionand so violates local injectivity (see Figure 2.12). If either feasibility condition isviolated we do not consider the candidate flip further. If, however, the flip is feasi-ble we evaluate energy decrease by checking the largest angle in the edge-adjacenttriangle pair both before and after the candidate flip is performed over all triangu-lations. If this largest angle decreases across meshes we flip the candidate edge toimprove local mesh quality.We apply edge split and edge collapse operations to bound edge lengths andangles to reasonable ranges while preserving our constraints. For each candidateedge, we consider an operation when it is 1. feasible with respect to our constraintsin Equation 2.8 2. out of relative bounds over all shape meshes, and 3. will locallyimprove the mesh. For our feasibility evaluation we check 1. if the operation leadsto an inversion and so violates local injectivity, 2. if the operation would deletea feature point (e.g., ruling out collapse of sharp corners), and 3. if the operationwould change the boundary topology (e.g., ruling out collapse of internal edgeswith both endpoints on the boundary). If any of these feasibility conditions isviolated we do not consider the candidate operation further. Feasible edge splitsare then applied when the edges relative length is too long in all meshes and theoperation would improve the worst edge length and local angles in the immediateneighborhood over all meshes. In turn, feasible edge collapses are applied when theedges relative length is too short in all meshes or if its opposite angles in incidenttriangles are too small. Please see Figure 2.12 for illustration and Section 4.7 foradditional evaluation of our Comesh Optimization solver.2.5 Interpolation EnergiesWith the well-formed multimeshes generated by our Comesh Optimization we cannow apply physically-motivated quantities as metrics for smooth, stable, and con-sistent interpolation over extreme changes in shape. We apply an elastic energy,WD(Xi,x), to measure deformation of shape x away from each example shape Xi,and introduce two new energies for interpolation, defined on shape boundaries.26With CohesiveElastic OnlyWith DynamicsFigure 2.14: Our cohesive energy integrates over cohesive zones (blue) that,together with sparse correspondences between shapes, resolve map-ping ambiguity to guide interpolation. Our cohesive energy and dy-namic blending enable adhesion and separation for topology changeas well as secondary physical effects in shape interpolation.We construct a cohesive energy, WC(Ci,x) to achieve topology changes for de-formed shape x across example shape cohesive edges Ci, and a non-overlap energy,WO(Bi,x), to avoid overlap of example shape boundaries Bi during variational in-terpolation, when desired. With these new energies shapes can now also separate,merge together, or even be prevented from overlapping during variational interpo-lation; see e.g., Figures 2.1, 2.14 and 2.16.To interpolate between shapes Xi, i ∈ [1,m], with topology changes and non-overlap, our full interpolation energy isWI(x,w) =m∑i=1wi[WC(Ci,x)+WO(Bi,x)+WD(Xi,x)], (2.9)with interpolation weights w=(w1, ...,wm)T ∈Rm+. Note that enforcing non-overlapis optional: when we wish to allow overlap in shape x, we simply omit the energyWO(·, ·) in Equation 2.9; see Figure 2.16. To find an interpolated shape x¯(w) wesolve Equation 2.3 with interpolation energy Equation 2.9.2.5.1 Cohesive EnergyWe introduce a new quadratic potential model, a cohesive energy, WC, that pullscuts and holes open and closed as we near example keyframes with differing topol-27ogy in our interpolation.Recall that cohesive edge sets, Ci, are collections of edge-edge pairs. Eachpair in shape i, given by C ji = {{va,vb},{vc,vd}} matches generally non-adjacentedges (Xai ,Xbi ) and (Xci ,Xdi ) for correspondence. They are created on the inte-rior and boundary of example shapes to indicate edge-pair correspondences wherean interpolated deformation will open as it moves towards example shapes whereedges are separated, and to merge together as the interpolation moves towards ex-ample shapes where cohesive edge-pairs are colocated. To do this we augment thestandard elastic measure of deformation in the shape average, WI , with an addi-tional cohesive energy, WC(Ci,x) for each example shape i, with cohesive edge set,Ci.Unlike an elastic potential, which integrates over the entire domain, our cohe-sive energy integrates solely along cohesive-zones. Thus each cohesive energyWC(Ci,x) integrates across the corresponding cohesive edge set Ci of exampleshape i. Figure 2.14 illustrates cohesive zones (blue) needed to separate individualshape (yellow and red; red and purple; yellow and purple) or to open a hole in theinterior (pink). Our cohesive potential measures the displacement difference witha quadratic form to enable efficient optimization with a local-global solver; seeSection 4.7.To construct WC we start by defining a quadratic measure of separation oncohesive edgesd(C ji ,x) = λ[xa− xcxb− xd]T [P 12 P12 P P][xa− xcxb− xd]. (2.10)Here λ ∈R+ enables control of cohesive stiffness and P is anR2×2 matrix that pro-vides control over topology change behavior to bias towards greater resistance toshear or stretch. Minimizing separation distance alone draws cohesive-zone edge-pairs together. To enable separation as deformation grows we then construct ourcohesive energy by clamping the quadratic, d, to an controlled maximum threshold280Figure 2.15: We measure displacment jump over cohesive edges using aclamped quadratic to model cohesive energy.τ; see Figure 2.15. Our cohesive energy is thenWC(Ci,x) = ∑C ji ∈Cimin( 13l ji d(Cji ,x), τ). (2.11)When WC is minimized in the shape average Equation 2.9, cohesive-zone edge-pairs are drawn towards each other when our distance measure is below τ , andseparate from one another as the threshold is exceeded. Here the multiplier 13 lji isobtained by piecewise-linear discretization of our energy (see our supplement fordetails) along cohesive edges C ji in shape Xi, withl ji = |Xai −Xbi |= |Xci −Xdi |. (2.12)Our shear and stretch control matrix P are then defined byP =α Id+(1−2α) xˆxˆT ∈ R2×2, withxˆ =xa+ xc− xb− xd‖xa+ xc− xb− xd‖ ,(2.13)and α ∈ [0,1] giving weighted control of the relative contributions of shearing to29stretching in the cohesive energy. Shearing effects becoming more significant aswe increase α . Here Id is the 2×2 identity matrix and xˆ is the edge-pair averageddirection defined on deformed shape x. Note that our cohesive energy then mod-els both separation and merging behaviors. Separated edge pairs that come closeenough to reduce the displacement field difference below thresholds τ are onceagain drawn together to merge.2.5.2 Non-overlap and Deformation EnergiesWhen interpolating it is often desirable to additionally avoid overlapping shapes.For this we need to ensure that boundary vertices Bi do not penetrate the interiorof example shape i in deformed shape, x. To do this we construct our non-overlapenergyWO(Bi,x) = ∑v j∈BiI(Bi,x(v j),x), (2.14)with the indicator functionI(Bi,x j,x) =∞, if x j ∈Ω(x,T,Bi),0, else. (2.15)Here Ω(x,T,Bi) gives the interior domain of our mesh at x, composed by the coverof deformed triangles in T , excluding boundary edges. The indicator function,I(Bi,x j,x), then correspondingly evaluates to an extreme penalty of ∞ if a bound-ary vertex x j is inside any triangle in the deformed configuration x. While stan-dard minimization techniques can be challenged by such nonsmooth objectives,within the local-global framework we employ, see Section 4.7, we can enforce ournonoverlap potential by direct projection in the local solve step [19].Deformation Energy Finally, To measure internal deformation, in this work, weuse the discrete, triangle-based as-rigid-as-possible (ARAP) energy [26, 79, 124] ,WD(Xi,x) = ∑tk∈Tκkak(Xi)‖Fk(Xi,x)−Rk(Xi,x)‖2F . (2.16)30With Non-overlapWithout Non-overlapFigure 2.16: Applying our non-overlap energy enables variational interpola-tion to preserve the concentric layout of this circles-in-circles example(bottom). Interpolation without our non-overlap energy ignores over-lap and so interpolation will then break the concentric layout (top).with Fk(Xi,x) retrieving the deformation gradient of triangle tk ∈ T with respectto deformed shape x and rest shape Xi. Here Rk is then the projection of Fk ontorotations and κk is the local material stiffness weighted per triangle, tk ∈ T , forfine-grained control of local resistance to deformation.2.6 Animating InbetweensWe have so far shown how to interpolate inbetweens from sets of arbitrarily drawnshapes. To build animations we can design sequences of inbetweens with an-imation curves w(s) ∈ Rm+, s ∈ R. Our animation function is then x¯(w(s)) =argminx WI(x,w(s)). Varying w(s) with s gives us control and exploration of timingand spacing, e.g., ease-in and ease-out, while explicitly treating s as a time variableenables the addition of physics-based dynamic effects to our inbetweens.2.6.1 Dynamic Effects for Variational InterpolationSubdividing the animation curve w(s), s ∈ [0,S] into sequential frames {x¯(w(0)),..., x¯(w(s)),.., x¯(w(S))}, we construct a discrete velocity for each frame with fi-nite differencing as x˙(s) = 1hs(x¯(w(s))− x¯(w(s−1))) where hs is the length of thedesired time interval between frames s and s−1.A new, time-stepped position for frame s due only to dynamics and thus without31any interpolation is thenxp(s) = x¯(w(s−1))+hsx˙(s−1)+h2s M−1 f (s). (2.17)Here M = ρI ∈R2n×2n is a lumped mass matrix with material density ρ while f (s)adds user specified forces, e.g., gravity, wind, damping and so forth, at time s.To blend dynamic effects into our interpolation we then construct an additionalenergy that biases towards a dynamics step for frame s, with a blending parameterξ ∈ R+,WT (x,ξ ) =ξ2h2s(x− xp(s))T M(x− xp(s)). (2.18)Minimized on its own, WT (x,ξ ) with respect to x this energy applies a forwardtime step of dynamic motion and ignores interpolation altogether.To inbetween with dynamic effects from frame s−1 to frame s, we then com-pose WT with the variational average to build a blending between variational inter-polation and dynamics effectsx¯(w(s),ξ (s)) = argminxWI(x,w(s))+WT (x,ξ (s)). (2.19)Here, the time varying blending parameter ξ (s) balances between interpolatedshapes and dynamics-driven deformation. When we set ξ (s) = 0, inbetweeningreturns us shapes x¯ that reproduce interpolated frames while, as ξ (s) grows, dy-namics pushes us farther from artist-designed shapes.While some artists may prefer to directly control the blending parameter ξ (s),we also provide a blending controller that automates the process of balancing be-tween dynamic effects and hitting keyframes in example shapes. Effectively, asthe difference between the dynamics and and keyframe shapes grow, our controlleradaptively drives down ξ (s). Our controller then allows greater physics-driven de-formation as it more coincides with drawn keyframe shapes. Artist-level controlis then simplified to a single guiding compliance term ε ∈ [0,1] that indicates howtightly inbetweening should follow keyframe shapes. See Appendix A for details.322.7 Experiments and EvaluationsIn this section we evaluate the performance of our Comesh Optimization solver;discuss and demonstrate our workflow for creating animations from examples; andpresent statistics and examples from animation synthesis with our methods. ForComeshing we examine our solver’s scaling behavior as we increase mesh resolu-tion, add more example shapes, and as we increase problem difficulty. For anima-tion applications, we examine our method by recreating the benchmark animationtasks from The 12 Basic Principles of Animation and on a range of challenginganimation examples with timing statistics.2.7.1 ImplementationAll evaluations and experiments are performed on a four-core Intel 3.50GHz CPUwith our C++ implementation of the system. Our implementation stores the inputmultimesh data in an augmented half-edge data structure that we use throughoutcomeshing and interpolation.For Comesh Optimization we adopt the fixed parameters in Algorithm 1 for allexamples. That is max outer itr and max inner itr are set to 500 and 100 re-spectively, while outer tol and inner tol are set to 10−6 and 10−7 respectively.For our mesh optimization and refinement we set angles and relative edge lengthbounds following Brochu et al. [21].For interpolation we apply local-global iterations [19, 79] to compute inbe-tweened animation frames. Our stopping criteria for the solver is set to terminatewhen the inf-norm of our relative residual is less than 10−3. The inbetweeningenergy we adopt is rigid motion invariant and we combine TRACKS solver [14]to resolve rigid frame during interpolation. In our examples we set our cohesivepotential stiffness λ = 4×103 while we vary τ . For dynamics we need not explic-itly set timesteps, h, material stiffnesses, κ , nor material density, ρ , as our dynamicblending model re-parameterizes these effects with the single user compliance termε; see our supplement for details and Table 2.3 for parameters used to finalize eachanimation example. Numerical entries for quadratic terms in the global matrix ofour interpolation solve, i.e., P in our cohesive energy Equation 2.11, are updatedin each global step. We then correspondingly update our linear system and numer-33ically factorize once per iteration. Finally, in each local step, along with standardARAP projections [19], wherever non-overlap energies are active intersections aresimply back-projected to the closest feasible point while, for the cohesive energy,edge-pair directions xˆ are updated.2.7.2 ComeshingMesh Resolution We first investigate the scaling of our Comesh Optimizationsolver’s performance as we increase problem size in terms of number of vertices.In Figure 2.17 (a) we test the scaling of our Comesh Optimization solver withsquare-to-hexagon examples increasing vertex count up to well over 300K verticesand see that our solver scales efficiently. Statistics for all examples in the plot arefully detailed in Appendix A. An example of which, at a mesh resolution of over17K vertices, is shown in Figure 2.20 left. Here our unoptimized solver convergesin under 4 minutes while achieving a low-distortion mapping with high-qualityelements.Number of Example Shapes We next investigate scaling behavior of our ComeshOptimization solver as we increase problem size in terms of the number of pro-Time (s)Number of Vertices5(a) mesh scale vs timingTime (s)Number of Shapes(b) mesh number vs timingFigure 2.17: Comesh Optimization statistics: our Comesh solver has efficientscaling as we increase both (a) mesh resolution and (b) number of ex-ample shapes used to build our interpolation.34vided example shapes In Figure 2.17 (b) we see that the Comesh solver also scalesefficiently as we increase the number of input shapes from 2 to 16 for our walkingcharacter example in Figure 2.20, right. As detailed in Appendix A, the solverremains convergent in under 2 minutes, while achieving a low-distortion mappingwith high-quality elements, see Figure 2.20 right.Problem Difficulty We then investigate the performance of our Comesh solver aswe increase problem complexity. In Figure 2.18 we illustrate a series of exampleswhere we increasingly thin out the spokes in the star shape example to increasecomplexity. As shown in the left column of Figure 2.18, this introduces increas-ingly more distorted and challenging Comesh Optimization problems. In the rightcolumn of Figure 2.18 we show the low-distortion mappings and high-quality ele-ments we compute with our Comesh solver, while in the top two rows of Figure 2.4we compare the improvement of our Comesh Optimized interpolations for hard 3and hard 6 examples over unoptimized LIM initializations. Finally, we observe thatsimilar input shapes generally require greater time to optimize; for timing statisticssee Appendix A. Currently, we conjecture that this may be due to similar inputshapes having greater ambiguity in conformal mapping and so requiring greatereffort to converge; we plan to investigate this further.Comesh Solver Statistics In Table 2.2 we detail the performance of our ComeshOptimization solver on our animation examples. Maximum iteration and timingwere 29 outer iterations with 18 seconds of compute time in our prototype code.In all example so far, with the exception of our alligator example, we observe thatComesh Optimization additionally reduces the number of triangles from the input,but not significantly, as its goal is refining topological connectivity of the initialmeshes.35Input OutputHard 1:Hard 2:26+4Hard 3:Hard 4:Hard 5:Hard 6:50 100 1505010015020025030050 100 15010020030040050050 100 1505010015020025030035050 100 15010020030040050060050 100 1505010015020025030050 100 15010020030040050050 100 1505010015020025030050 100 15010020030040050 100 1505010015020025030050 100 1505010015020025030035050 100 1505010015020025030050 100 15050100150200250Angle Distribution Angle DistributionFigure 2.18: Comesh Optimization examples over increasing problem com-plexity. Here our Comesh solver consistently produces low-distortionmappings and high-quality elements as we increase variation betweeninput shapes. We visualize input/output map distortions and angle dis-tributions here and compare the improvement of our resulting ComeshOptimized interpolations for hard 3 and hard 6 examples over unopti-mized LIM initializations in Figure 2.4.36Table 2.2: We list the comesh optimization statistics for our inbetweening examples: number of example shapes, in-put/output mesh information like number of vertex, triangle and cohesive edge; inner iteration number shows thetotal loops of geometrical and topological optimization while outer iteration number shows the number of outerloops; flip, collapse and split timing show the amount of time spent on respective topological operation while totaltiming shows the overall cost of the whole algorithm.Shape(#) Ver(#/#) Tri(#/#) CoE(#/#) Inner(#) Outer(#) Flip(s) Collapse(s) Split(s) Total(s)aligator 4 297/358 485/593 174/162 1051 19 1.722 1.615 0.857 7.589cell splitting 3 573/570 978/972 142/142 572 15 1.009 0.547 0.487 4.078chinese character 2 277/277 382/379 10/10 36 4 0.023 0.012 0.012 0.119circle in circle 2 201/201 241/241 126/126 70 3 0.024 0.012 0.013 0.137continental drift 4 271/269 362/357 104/104 2565 29 8.841 5.539 1.067 18.322growing flower 4 163/159 213/205 24/24 54 5 0.023 0.018 0.028 0.104flying bird 3 226/226 265/265 0/0 60 7 0.028 0.018 0.016 0.141stylized walk 16 132/132 161/161 0/0 243 6 0.311 0.166 0.166 0.824walk star 8 89/89 93/93 24/24 243 7 0.085 0.049 0.051 0.28437Table 2.3: We list the animation statistics and parameter settings for the inbetweening examples, including maximaliteration number per frame, maximal time spent on a single local-global step, compliance parameter ε , cohesivezone parameter τ , model control α , damping paramter µ as well as material type (homogeneous or inhomogeneous)and with or without animation curve tuning.Iteration(#) Local Step(s) Global Step(ms) ε τ α µ Uniform Material Animation Curvealigator 47 0.194 2.113 0.02 1e5 0.7 0.5 3 3cell splitting 173 0.202 6.144 1.7 2e6 0.2 1.0 3 7chinese character 38 0.133 1.307 0 2e5 0.5 - 3 7circle in circle 52 0.007 2.106 1.5 4e5 0.5 1.8 3 7continental drift 79 0.085 1.817 0 2e5 0.5 - 3 7growing flower 53 0.089 0.592 0 2e5 0.5 - 3 7flying bird 43 0.052 0.502 1 - - 0.5 7 3stylized walk 40 0.038 0.441 0 - - - 3 3swirl 406 0.04 0.273 0 - - - 3 7walk star 36 0.027 0.286 0 6e5 0.2 - 3 3382.7.3 AnimationControl Our system exposes a small number of parameters for interactive explo-ration of animations. These terms are listed in Table 2.3 and are summarized here.The blending compliance parameter ε , indicates the desired balance between dy-namic effects and example shape conformance; it is used to control how tightlyinterpolated dynamics match keyframe shapes; smaller values give interpolationsthat matches the keyframes better while larger value bias towards more exaggeratedynamics. The cohesive threshold, τ , controls the effective stickiness of separationand merging during interpolation. Our mode control parameter, α , then determineswhether stretch or shear dominates the topology change. The damping force scal-ing term, µ , is enabled when dynamics are applied and helps guide how excitedsecondary dynamics remain over time. Finally, we enable the ability to paint mul-tiple material stiffnesses, κ , with a brush interface to control the degree to whichshapes resist deformation.Range Our method generates interpolated animations for a range of challenginganimation problems consisting of extreme deformation, topology changes and dy-namics. Some of these examples are illustrated in Figure 2.1. Our system enablesthe generation of a wide variety of animations from sparse input drawings. As abaseline test we reproduce a set of benchmark examples from the 12 Basic Princi-ples of Animation, as shown in Appendix A. For this task we use a small numberof input drawings; always less than five. We additionally duplicate a professional,Figure 2.19: Trumperfly: interpolation from a sparse set of production framesenables experimentation and variation.39Example ExtendExample Scale26+4Input:Output:50 100 150100020003000400050 100 15010002000300040005000600050 100 150500100015002000250050 100 15020040060080010001200Angle Distribution Angle DistributionFigure 2.20: Our Comesh Optimization algorithm scales well to large sizemeshes (left) and large number of inputs (right). In both cases, ourmethod efficiently generates meshes with high-quality elements andmaps.hand-drawn animation that includes extreme transitions with topology change andextreme deformations; see Figure 2.19 for sample interpolated frames. We usejust a sparse subset of the frames and then demonstrate freeform experimentationon the resulting interpolation to interactively explore animation variations withouthaving to draw new frames.Statistics Our interpolation provides interactive feedback so that artists can quicklyiterate in creating final animations. Performance statistics for interpolating our an-imation examples (both with and without dynamics) are summarized in Table 2.3with a maximum of less than 10ms per local-global iteration across all examples inour unoptimized code.We have presented an interpolation system that allows artists to interactivelyanimate and interpolate inbetweens from arbitrary key drawing shapes. We demon-strate that this framework enables exploration of a wide range of 2D animationsfrom sparse drawing inputs. Our method provides artists the flexibility to specifydesired correspondences and topology changes, and to create and edit inbetweeningfor animation. To achieve artifact-free interpolation, we introduce a Comesh Op-timization method and solver that builds consistent meshes across drawings whilerespecting correspondences and cuts. Our optimization algorithm improves both40the energy discretization on each mesh, by enhancing element quality, as wellas the mapping distortion across meshes. Given our Comesh Optimization, wethen introduce two boundary energies, a cohesive energy and a non-overlap en-ergy. Finally, we demonstrate how to add dynamic secondary motion effects to ourkeyframe interpolation.Looking ahead we believe that a complete system for drawing planar animationshould also support layered drawings. A few of our examples use two layers (nearand far) and rely on interpolation to preserve coherence. More elaborate drawingarrangements will likely be enhanced by multilayer support: techniques to extractlayers from drawings and techniques to attach them together explicitly. ComeshOptimization could also be further improved with a smooth spline-defined bound-ary. As is, a feature point, such as a sharp corner, on the boundary will not moveduring optimization. Likewise, we expect that the numerical solver for our ComeshOptimization could be enhanced with improved convergence.Our approach could be extended to 3D. Our main contributions, Comesh Op-timization and our Multimesh Interpolation with boundary energies appear to gen-eralize quite easily to 3D. However, we also depend on topology optimization andpiecewise linearly injective mappings, which remain challenging to extend reliablyand efficiently in 3D.We note that in this work we do not address the interesting and important ques-tion of how to meaningfully automate correspondences for artists. Instead, wesought to expand the expressive range of interpolation so that artists can accom-plish more with less, such as complex animation of cell splitting with only threekeyframes. In doing so, we also recognized that correspondence plays a key se-mantic role and designed a solution that could work for correspondences discov-ered by either interactive or automatic workflows. As detailed here and illustratedin our inbetweening examples, our current method uses interactively defined cor-respondences from artists. Optimal transport based methods are one interestingavenue of future exploration for automation. Currently, however, as we illustratein Figure 2.2, such methods generally choose one path and not always the desiredone. Thus we see an investigation of automated and semi-automated methods forproviding correspondences remains an avenue of exciting future work here.Finally, the projected gradient descent strategy adopted in the inner loop of41our Comesh Optimization framework provides satisfying results but suffers fromslow convergence. In the following two chapters, we will investigate approachesto speed up the mapping solver. Next, we first turn to nonlinear optimization tech-niques and propose a new solver with state-of-the-art performance.42Chapter 3Blended Cured Quasi-Newton forDistortion Optimization3.1 IntroductionMany fundamental physical and geometric modeling tasks reduce to minimizingnonlinear measures of distortion over meshes. Simulating elastic bodies, defor-mation, parametrization, shape interpolation, deformable inverse kinematics, andanimation all require robust, efficient, and easily automated distortion optimiza-tion. By robust we mean the algorithm should solve every reasonable problemto any accuracy given commensurate time, and only report success when the accu-racy has truly been achieved. By efficient we mean rapid convergence in wall-clocktime, even if that may mean more (but cheaper) iterations. By automated we meanthe user needn’t adjust algorithm parameters or tolerances at all to get good resultswhen going between different problems. With these three attributes, a distortionoptimization algorithm can be reliably used in production software.We propose a new algorithm, Blended Cured Quasi-Newton (BCQN), withthree core contributions based on analysis of where prior methods faced difficulties:• an adaptively blended quadratic energy proxy for distortion energies thatiteratively combines the Sobolev gradient and a quasi-Newton secant ap-proximation; this allows low cost per iterate with second-order acceleration430 1000 2000 3000 4000 5000 6000 70000102030405060BCQNAQPPN0 1000 2000 3000 4000 5000 6000 70000102030405060201/207/3BCQN Converged0 1000 2000 3000 4000 5000 6000 700001020304050600 1000 5000 7000Rest ShapeInitializerTime (seconds)9/13/0 converged/550/9 converged/685/1270/77/1 137/138/2170/172/2160/162/2Figure 3.1: Twisting. A stress-test 3D deformation problem. Left: weinitialize a 1.5M element tetrahedral mesh bar with a straight restshape into a tightly twisted coil, constraining both ends to stay fixed.Right: minimizing the ISO deformation energy to find a constrainedequilibrium with (top to bottom) Projected Newton (PN), AcceleratedQuadratic Proxy (AQP) and our Blended Cured Quasi-Newton (BCQN)method, we show intermediate shapes at reported wall-clock time (sec-onds) and iteration counts at those times (BCQN/AQP/PN). AQP, muchslower than BCQN, requires many more iterations to converge whilePN, despite requiring fewer iterations, is well over 25X slower due toper-iteration costs dominated by factorization.while avoiding secant artifacts where the Laplacian is more robust;• a barrier-aware filter on search directions, that gains larger step sizes andso improved convergence progress in line search for the common case ofiterates where individual elements degenerate towards collapse; and• a characteristic gradient norm convergence criterion, immune to termi-nating prematurely due to algorithm stagnation, that is consistent acrossmesh sizes, scales, and choice of energy so that per-problem adjustment isunnecessary.Over a wide range of test cases we show that BCQN achieves our goals for produc-tion software: BCQN makes the solution of some previously intractable problemspractical, offers up to an order of magnitude speed-up in other cases (see Fig-ure 3.1) and, in all cases investigated so far, robustly converges while improvingon or closely matching the performance of the best-in-class optimizers available.443.2 Problem Statement and OverviewDistortion optimization seeks local minimizersx∗ = argminx∈RdnE(x), (3.1)for n vertex locations in d-dimensional space stored in vector x, where the energyE(x) is a measure of distortion, and x is subject to constraints.1 The energy isexpressed as a sum over elements t in a triangulation T (triangles or tetrahedradepending on dimension),E(x) = ∑t∈TvtW(Ft(x)), (3.2)where vt > 0 is the area or volume of the rest shape of element t, W is an energydensity function taking the deformation gradient as its argument, and Ft computesthe deformation gradient for element t. This problem may be given as is, or maybe the result of a discretization of a continuum problem with, for example, linearfinite elements.3.2.1 Iterative Solvers for Nonlinear MinimizationSolution methods for distortion optimization generally apply an algorithmic strat-egy of iterated approximation and stepping [15], built from three primary ingredi-ents: an energy approximation, a line search, and a termination criteria.2Energy ApproximationAt the current iterate xi we form a local quadratic approximation of the energy, orproxy:Ei(x) = E(xi)+(x− xi)T∇E(xi)+ 12(x− xi)T Hi(x− xi) (3.3)1For simplicity of presentation we restrict our attention here to constraining a subset of vertexpositions to given values, i.e. Dirichlet boundary conditions.2Alternatively, trust-region methods are available, though not considered in the current work noras popular within the field.45where Hi is a symmetric matrix. Near the solution, if Hi accurately approximatesthe Hessian we can achieve fast convergence optimizing this proxy, but it is alsocritical that it be stable — symmetric positive definite (SPD) — to ensure the proxyoptimization is well-posed everywhere; we also want Hi to be cheap to solve with,preferring sparser matrices and ideally not having to refactor at each iteration.Line SearchQuadratic models allow us to iteratively apply linear solvers to find stationarypoints x∗i = argminx Ei(x) of the local energy approximation. A steppi = xi∗− xi =−H−1i ∇E(xi) (3.4)towards this stationary point then forms a direction for probable energy descent.However, quadratic models are only locally accurate for nonlinear energies in gen-eral, thus line-search is used to find an improved length αi > 0 along pi to updateto a new iteratexi+1← xi+αi pi, (3.5)that provides adequate decrease in nonlinear energy E. Of particular concern forthe distortion problems we face are energies that blow up to infinity for degenerate(flattened) elements: in a given step, the elements where this may come close tohappening rapidly depart from the proxy, and the step size αi may have to be verysmall indeed, see Figure 3.2, impeding progress globally.TerminationIteration continues until we are able to stop with a “good enough” solution – butthis requires a precise computational definition. Typically we monitor some quan-tity that approaches zero if and only if the iterates are approaching a stationarypoint. The standard in unconstrained optimization is to check the norm of the gra-dient of the energy, which is zero only at a stationary point and otherwise positive;46Descent Direction Step After Line Search PNLine SearchBCQNAQPEmbedding SurfaceFigure 3.2: Line-search filtering. Barrier terms in nonconvex energies (herewe use ISO) of the form 1/g(σ) can severely restrict step sizes in linesearches even when expensive, high-quality methods such as Newton-type methods are applied. Left column: descent-direction vector fields,per vertex, in a descent step generated by BCQN, PN and AQP withpotential blocking triangles rendered in red. Right column, bottomrows: after line-search, close to collapsing elements have restricted theglobal step size for AQP and PN that block progress. Right column,top row: BCQN’s barrier-aware line-search filtering cures the descentdirection enabling progress with significant descent.47however, the raw gradient norm depends on the mesh size, scaling, and choice ofenergy, which makes finding an appropriate tolerance to compare against highlyproblem-dependent and difficult to automate.3.3 Related Work3.3.1 Energies and ApplicationsA wide range of physical simulation and geometry processing computations arecast as variational tasks to minimize measures of distortion over domains.To simulate elastic solids with large deformations we typically minimize hyper-elastic potentials formed by integrating strain-energy densities over the body. Thesematerial models date back to Mooney [86] and Rivlin [110]. Their Mooney-Rivlin and Neo-Hookean materials, and many subsequent hyperelastic materials,e.g. St. Venant-Kirchoff, Ogden, Fung [18], are constructed from empirical ob-servation and analysis of deforming real-world materials. Unfortunately, all but afew such energy densities are nonconvex with respect to strain.3 This makes theirminimization highly challenging. Constants in these models are determined by ex-periment for scientific computing applications [97], or alternately are directly setby users in other cases [150], e.g., to meet artistic needs.In geometry processing a diverse range of energies have been proposed to min-imize various mapping distortions, generally focused on minimizing either mea-sures of isometric [1, 25, 78, 120, 123] or conformal [12, 37, 53, 73, 87, 140]distortion. While some of these energies do not prohibit inversion [25, 37, 73, 123]many others have been explicitly constructed with nonconvex terms that guaranteepreservation of local injectivity [1, 53, 120]. Other authors have also added con-straints to strictly bound distortion [17, 76] for example, but we restrict attention tounconstrained minimization — and note constrained optimization often relies onunconstrained algorithms as an inner kernel.Our goal here is to provide a tool to minimize arbitrary energy density func-3We observe this is distinct from the global energy being nonconvex w.r.t. vertex positions:ARAP [123] is an example of a convex energy density that nonetheless can lead to a nonconvexglobal energy.48tions as-is. We take as input energy functions provided by the user, irrespectiveof whether these energies are custom-constructed for geometry tasks, physical en-ergies extracted from experiment, or energies hand-crafted by artists. Our workfocuses on the better optimization of the important nonconvex energies whoseminimization remains the primary challenging bottleneck in many modern geom-etry and simulation pipelines. In the following sections, to evaluate and comparealgorithms, we consider a range of challenging nonconvex deformation energiescurrently critical in physical simulation and geometry processing: Mooney-Rivlin(MR) [20], Neo-Hookean (NH) [20], Symmetric Dirichlet (ISO) [120], ConformalDistortion (CONF) [1], and Most-Isometric Parameterizations (MIPS) [53].3.3.2 Energy ApproximationsBroadly, existing models for the local energy approximation in Equation 3.3 fallinto four rough categories that vary in the construction of the proxy4 matrix Hi.Newton-type methods exploit expensive second-order derivative information; first-order methods use only first derivatives and apply lightweight fixed proxies; quasi-Newton methods iteratively update proxies to approximate second derivatives usingjust differences in gradients; Geometric Approximation methods use more domainknowledge to directly construct proxies which relate to key aspects of the energy,resembling Newton-type methods but not necessarily taking second derivatives.Newton-type methods generally can achieve the most rapid convergence butrequire the costly assembly, factorization and backsolve of new linear systems perstep. At each iterate Newton’s method uses the energy Hessian, ∇2E(xi), to forma proxy matrix. This works well for convex energies but requires modification fornonconvex energies [93] to ensure that the proxy is at least positive semi-definite(PSD). Composite Majorization (CM), a tight convex majorizer, was recently pro-posed as an analytic PSD approximation of the Hessian [115]. The CM proxy isefficient to assemble but is limited to two-dimensional problems and just a trio ofenergies: ISO, NH and symmetric ARAP. More general-purpose solutions includeadding small multiples of the identity, and projection of the Hessian to the PSD4Names and notations for Hi vary across the literature depending on method and application. Forconsistency, here, across all methods we will refer to Hi as the proxy matrix — inclusive of caseswhere it is the actual Hessian or direct modification thereof.49cone but these generally damp convergence too much [80, 93, 115]. More effec-tive is the Projected Newton (PN) method that projects per-element Hessians to thePSD cone prior to assembly [134]; also of interest is Chen and Weber’s analyticHessian projection in a reduced basis setting [27]. Both CM and PN generallyconverge rapidly in the nonconvex setting with CM often outperforming PN in thesubset of 2D cases where CM can be applied [115], while PN is more generalpurpose for 3D and 2D problems. Both PN and CM have identical per-elementstencils and so identical proxy structures. Despite low iteration counts they bothscale prohibitively due to per-iteration cost and storage as we attempt increasinglylarge optimization problems.First-order methods build descent steps by preconditioning the gradient with afixed proxy matrix. These proxies are generally inexpensive to solve and sparseso that cost and storage remain tractable as we scale to larger systems. How-ever, they often suffer from slower convergence as we lack higher-order informa-tion. Direct gradient descent, Hi← Id, and Jacobi-preconditioned gradient descent,Hi← diag(∇2E(xi))offer attractive opportunities for parallelization [41, 138] butsuffer from especially slow convergence due to poor scaling. The Laplacian ma-trix, L, forms an excellent preconditioner, that both smooths and rescales the gra-dient [70, 82, 92]. Unlike the Hessian, the Laplacian is a constant PSD proxy thatcan be pre-factorized once and backsolved separately per-coordinate. Iterating de-scent with Hi← L, is the Sobolev-preconditioned gradient descent (SGD) method.SGD was first introduced, to our knowledge, by Neuberger [92], but has since beenrediscovered in graphics as the local-global method for minimizing ARAP [123].As noted by Kovalsky et al. [70] local-global for ARAP is exactly SGD. More re-cently Kovalsky et al. [70] developed the highly effective Accelerated QuadraticProxy (AQP) method by adding a Nesterov-like acceleration [91] step to SGD.This improves AQP’s convergence over SGD. However, as this acceleration is ap-plied after line search, steps do not guarantee energy decrease and can even containcollapsed or inverted elements — preventing further progress; see e.g. Figure 3.4.More generally, the Laplacian is constant and so ignores valuable local curvatureinformation — we see this issue in a number of examples in Section 3.8 whereAQP stagnates and is prohibitively slow to converge. Curvature can make the crit-ical difference to enable progress.50Quasi-Newton methods lie in between these two extremes. They successively,per descent iterate, update approximations of the system Hessian using a varietyof strategies. Quasi-Newton methods employing sequential gradients to updatesproxies, i.e. L-BFGS and variants, have traditionally been highly successful inscaling up to large systems [15]. Their updates can be performed in a computeand memory efficient manner and can guarantee the proxy is SPD even where theexact Hessian is not. While not fully second-order, they achieve superlinear con-vergence, regaining much of the advantage of Newton-type methods. L-BFGSconvergence can be improved with the choice of initializer. Initializing with the di-agonal of the Hessian [93], application-specific structure [59] or even the Laplacian[80] can help. However, so far, for distortion optimization problems, L-BFGS hasconsistently and surprisingly failed to perform competitively irrespective of choiceof initializer [70, 108]. Nocedal and Wright point out that the secant approxima-tion can implicitly create a dense proxy, unlike the sparse true Hessian, directlyand incorrectly coupling distant vertices. This is visible as swelling artifacts forintermediate iterations of the L-BFGS methods in Figure 3.5.Geometric Approximation methods for distortion optimization have also beendeveloped recently: SLIM [108] and the AKAP preconditioner [32]. SLIM ex-tends the local-global strategy to a wide range of distortion energies while AKAPapplies an approximate Killing vector field operator as the proxy matrix. Bothrequire re-assembly and factorization of their proxies for each iterate. SLIM andAKAP convergence are generally well improved over SGD and AQP [32, 108].However, they do not match the convergence quality of the second-order, Newton-type methods, CM and PN [115]. SLIM falls well short of both CM and PN [115].AKAP is more competitive than SLIM but remains generally slower to convergethan PN in our testing, and is much slower than CM. At the same time SLIM andAKAP stencils, and so their fill-in, match CM’s and PN’s; see Figure 3.8. SLIMand AKAP thus require the same per-iteration compute cost and storage for linearsolutions as PN and CM without the same degree of convergence benefit [115].In summary, for smaller systems Newton-type methods have been, till now,our likely best choice for distortion optimization, while as we scale we have in-evitably needed to move to first-order methods to remain tractable, while acceptingreduced convergence rates and even the possibility of nonconvergence altogether.51We develop a new quasi-Newton method, BCQN, that adaptively blends gradientinformation with the matrix Laplacian at each iterate to regain improved and robustconvergence with efficient per-iterate storage and computation across scales whileavoiding the current pitfalls of L-BFGS methods.3.3.3 Line SearchOnce we have applied the effort to compute a search direction we would like tomaximize its effectiveness by taking as large a step along it as possible. Becausethe energies we treat are nonlinear, too large a step size will actually make thingsworse by accidentally increasing energy. A wide range of line-search methodsare thus employed that search along the step direction for sufficient decrease [93].However, when we seek to minimize nonconvex energies on meshes the situa-tion is more challenging. Most (although not all) popular and important nonlin-ear energies, both in geometry processing and physics, are composed by the sumof rational fractions of singular values, σ , of the deformation gradient W (F) =W (σ) = f (σ)/g(σ) where the denominator g(σ)→ 0 as σi→ 0,∀i ∈ [1,d]. No-tice that these 1/g(σ) barrier functions block element inversion. Irrespective oftheir source, these blocking nonconvex energies rapidly increase energy along anysearch direction that would collapse elements. To prevent this (and likewise thepossibility of getting stuck in an inverted state) search directions are capped to pre-vent inversion of every element in the mesh. This is codified by Smith and Scha-effer’s [120] line-search filter, applied before traditional line search, that computesthe maximal step size that guarantees no inversions anywhere.Unfortunately, this has some serious consequences for progress. Notice that ifeven a single element is close to inversion this can amputate the full descent stepso much that almost no progress can be made at all; see Figure 3.2, bottom andmiddle rows. This, in many senses seems unfair as we should expect to be able tomake progress in other regions where elements may be both far from inversion andyet also far from optimality. To address these barrier issues we develop an efficientbarrier-aware filter that allows us to avoid blocking contributions from individualelements close to collapse while still taking large steps elsewhere in the mesh, seeFigure 3.2, top row.523.3.4 TerminationNaturally we want to take as few iterates as possible while being sure that whenwe stop, we have arrived at an accurate solution according to some easily specifiedtolerance. The gold-standard in optimization is to iterate until the gradient is small‖∇E‖ < ε , for a specified tolerance ε > 0. This is robust as ∇E is zero only atstationary points, and with a bound on Hessian conditioning near the solution caneven provide an estimate on the distance of x to the solution.Coarse meshFine mesh1e-41e-3ToleranceFigure 3.3: Standard termination measures, e.g. the vertex-scaled gradientnorm above, are inconsistent across mesh, energy and scale changes.However, an appropriate value of ε for a given application is highly dependenton the mesh, its dimensions, degree of refinement, energy, etc. A common engi-neering rule of thumb to deal with refinement consistency is to instead divide theL2-norm of ∇E by the number of mesh vertices. However, as we see in Figure 3.3,this normalization does not help significantly, for example here across changes inmesh resolution for the 2D swirl test; see Section 3.8.2 for more experiments.To avoid problem dependence, recent distortion optimization codes generallyeither take a fixed (small) number of iterations [108] or iterate until an absolute orrelative error in energy |Ei+1−Ei| and/or position ‖xi+1− xi‖ are small [70, 115].However, experiments underscore there is not yet any method which always con-verges satisfactorily in the same fixed number of iterations across varying boundaryconditions, shape difficulty, mesh resolution, and choice of energy. Measuring the53change in energy or position, absolutely or in relative terms, unfortunately cannotdistinguish between an algorithm converging and simply stagnating in its progressfar from the solution; again, there is not yet any method which can provably guar-antee any degree of progress at every iterate before true convergence. Figure 3.4illustrates on the swirl example how a reference AQP implementation declares con-vergence well before it reaches a satisfactory solution, when early on it hits a diffi-cult configuration where it makes little local progress.To provide reassuring termination criteria in practice and to enable fair compar-isons of current and future distortion optimization problems we develop a gradient-based stopping criterion that remains consistent for optimization problems even aswe vary scale, mesh resolution and energy type. This allows us, and future users,to set a default convergence tolerance in our solver once and leave it unchanged,independent of scale, mesh and energy. This likewise enables us to compare al-gorithms without the false positives given by non-converged algorithms that havehalted due to lack of progress.BCQNAQPFigure 3.4: In the 2D swirl example, BCQN with our reliable terminationcriterion (right) only stops once it has actually reached a satisfactorysolution. The reference AQP implementation (left) erroneously declaressuccess early on when it finds two iterates have barely changed, but thisis due only to hitting a difficult configuration where AQP struggles tomake progress.3.4 Blended Quasi-NewtonIn this section we construct a new quadratic energy proxy which adaptively blendsthe Sobolev gradient with L-BFGS-style updates to capture curvature information,avoiding the troubles previous quasi-Newton methods have encountered in dis-tortion optimization. Apart from the aforementioned issue of a dense proxy in-54correctly coupling distant vertices in L-BFGS, we also find that the gradients fornon-convex energies with barriers can have highly disparate scales, causing furthertrouble for L-BFGS. The much smoother Sobolev gradient diffuses large entriesfrom highly distorted elements to the neighborhood, giving a much better scaling.The Laplacian also provides essentially the correct structure for the proxy, onlydirectly coupling neighboring elements in the mesh, and is well-behaved initiallywhen far from the solution, thus we seek to stay close to the Sobolev gradient, asmuch as possible, while still capturing valuable curvature information from gradi-ent history.The standard (L-)BFGS approach exploits the secant approximation from thedifference in successive gradients, yi =∇E(xi+1)−∇E(xi) compared to the differ-ence in positions si = xi+1− xi,∇2E(xi+1)si ' yi⇒ ∇2E(xi+1)−1yi ' si,(3.6)to update an inverse proxy matrix Di =H−1i (approximating ∇2E−1 in some sense)so that Di+1yi = si. The BFGS quasi-Newton update is genericallyQNi(z,D) =Vi(z)T DVi(z)+sisTisTi z, Vi(z) = I− zsTisTi z. (3.7)We can understand this as using a projection matrix Vi to annihilate the old D’saction on z, then adding a positive semi-definite symmetric rank-one matrix toenforce QNi(z,D)z = si. Classic BFGS uses Di+1 = QNi(yi,Di), whereas L-BFGSusesDi+1 = QNi(yi, D˜i), (3.8)where D˜i has the oldest QN update removed, and crucially represents each D asa product of linear operators, rather than an explicit full matrix. Only the last m{s,y} vector pairs (we use m = 5) along with the initial D1 (we use the inverseLaplacian, storing only its Cholesky factor) are stored; application of D is then justa few vector dot-products and updates along with backsolves for the Laplacian.5501002003004004.0.55.0.56.0.57.0.5SGDL-BFGSSL-BFGSAQPBL-LBFGSBCQN100 200 300 4004.04.55.05.56.06.57.07.5BCQN Blended AQPIterate 0SL-BFGSIterate 240L-BFGS SGDlendedIterationsFigure 3.5: A 2D shearing deformation stress test with MIPS energy, com-paring methods by plotting iteration vs. energy. Both L-BFGS as wellas Laplacian initialized L-BFGS (SL-BFGS) have slow convergence aspreviously reported [70] — especially when compared to SGD and AQPwhich use just the Laplacian. At iteration 240 the visualized deforma-tions show both L-BFGS-based methods suffering from swelling due toinaccurate coupling of distant elements. Applying our blending modelalone (Blended) is highly effective, while our full BCQN method givesthe best results overall.3.4.1 Greedy Laplacian BlendingExperiments show that far from the solution, the Laplacian is often (althoughnot always) a much more effective proxy than the L-BFGS secant version: seeAQP/SGD vs. L-BFGS in Figure 3.5. In particular, the difference in energies ymay introduce spurious coupling or have badly scaled entries near distorted trian-gles. In this case if the energy were based on the Laplacian itself (the Dirichletenergy), the difference in gradients would be the better behaved Ls. This motivatestrying to update with Ls instead of y,Di+1 = QNi(Lsi, D˜i). (3.9)56This would keep us consistent with Sobolev preconditioning, which is often veryeffective in initial iterations. However, to achieve the superlinear convergence L-BFGS offers, near solutions we wish to come closer to satisfying the secant equa-tion, and so aim to move towards using y instead. More generally we seek tobalance between gradient information and the action of the Laplacian.We thus propose an adaptive blending strategy that updates the proxy per iterateusingzi = (1−βi)yi+βiLsi (3.10)in QN(zi, D˜i), with a blending parameter βi ∈ [0,1] constructed below. Towardsour goals we first derive here a greedy blending that is “as Laplacian as possible”to approximate the energy gradient. Then, below in Section 3.4.2, we present ourfinal, empirically derived variant that we have experimentally determined to so farwork best. As we will soon see this blend is highly effective; however, further workremains to find and confirm optimal blending terms.We begin by greedily selecting a βi to scale Lsi to be as close to yi as possible,βi = argminβ∈[0,1]‖yi−βLsi‖2, (3.11)in other words using the projection of yi onto Lsi. This comes as close as possibleto satisfying the secant equation with Lsi and then makes up the remainder with yi.Solving Equation 3.11 givesβi = proj[0,1](yiT Lsi‖Lsi‖2). (3.12)Observe that when Ls is roughly aligned with the gradient jump y , but y is as largeor larger, β grows and Laplacian smoothing increases — as we might hope forinitially when far from the solution, where the Sobolev gradient is most effective.When the energy Hessian diverges strongly from the Laplacian approximation, per-haps when the cross-terms between coordinates missing from the scalar Laplacianare important, then β will decrease, so that contributions from yi again grow. Fi-nally, as the gradient magnitudes decreases close to the solution, β will similarlydecay, ideally regaining the superlinear convergence of L-BFGS near local minima.573.4.2 Blended Quasi-NewtonWith the blending projection Equation 3.12 in place we experimented with a rangeof rescalings in hopes of further improving efficiency and robustness. After ex-tensive testing we have so far found the following scaling to offer the best perfor-mance:βi = proj[0,1](normest(L)yiT LsiB(T )),with B(T ) =(∑t∈Tvt) 2(d−1)d.(3.13)Here normest(L) is an efficient estimate of the matrix 2-norm using power itera-tion [51], and B(T ) is a constant normalizing term with appropriate dimensionsand so no longer has the same potential concern for sensitivity in the denominatorwhen Ls is small but s is not. Both terms are computed just once before iterationsbegin and reused throughout.As mentioned, we initialize the inverse proxy with D1 = L−1, thus startingwith Laplacian preconditioning. With line search satisfying Wolfe conditions ourproxy remains SPD across all steps [93] and so delivers descent directions. Eachstep jointly updates Di using the standard two-loop recursion and finds the nextdescent direction si = −Di∇E(xi). Figure 3.5 illustrates the gains possible fromblended quasi-Newton compared to both standard L-BFGS and Sobolev gradientalgorithms, while then applying our barrier-aware filter, derived in the next sectiongives best results with our full BCQN algorithm.3.5 Barrier-Aware Line Search FilteringAs discussed in Section 3.3.3 and shown in Figure 3.2, barrier factors 1/g(σ) innonconvex energies typically dominate step size in line search. Even a single ele-ment that is brought close to collapse by the descent direction, pi, can restrict theline search step size severely. The computed step size αi then scales pi globally sothat all elements, not just those that are going to collapse along pi, are preventedfrom making progress. To avoid this, a natural strategy suggests itself: when thedescent direction would cause elements to degenerate towards collapse along the58Rest ShapeInitial Shape (a) Filtered Direction (b) Filtered Gradient (c) BCQNFigure 3.6: Direct filtering does not work. Zeroing out inverting compo-nents of descent directions or gradients makes the search direction in-consistent with the objective and so prevents convergence, leading totermination at poor solutions (a) and (b). Left: we initialize a 2D sheardeformation, constraining the top of a bar to slide rightwards. Middle:direct filtering of the descent direction (a) and the gradient (b) allowlarge descent steps forward unblocked from the contributions of close-to-collapsed elements. However, this results in termination at shapesthat that do not satisfy optimality of the original minimization. Right:compare to a locally optimal solution for this problem (c) obtained withBCQN.full step, rather than simply truncating line search as in Smith and Schaefer [120],we filter collapsing contributions from the search direction prior to line search. Wecall this strategy barrier-aware line search filtering.3.5.1 Curing Line SearchFigure 3.6 illustrates how simple, direct filters that zero out contributions fromnearly-inverted elements in either the search direction Figure 3.6a or the gradientbefore Laplacian smoothing Figure 3.6b fail. We must be able to make progress innearly-inverted elements when the search direction can help, or there is no hope forreaching the actual solution; simple zeroing fails to converge, which is no surpriseas it in essence is arbitrarily manipulating the target energy, changing the problembeing solved. We instead seek to augment the original optimization problem ina way that does not change the solution, but gives us a tool to safely deal withproblem elements so that the resulting search direction pi does not cause them tocollapse or invert, ideally with a small fixed cost per iteration.593.5.2 One-sided Barriers in Distortion OptimizationElement t ∈ T is inverted at positions x precisely when the orientation functionat(x) = det(Ft(x)) is negative. Concatenating over T , the global vector-valuedfunction for element orientations is thena(·) = (a1(·), ...,am(·))T . (3.14)As long as a(x) > 0, no element is collapsed or inverted, and the energy remainsfinite. Note, however, many energies are also finite for inverted elements at(x)< 0,only blowing up at collapse at(x) = 0, so technically there may exist local minimawhere ∇E(x∗) = 0 yet some elements are inverted. Generally, practitioners wish torule these potential solutions out. However, this is done with two implicit, but sofar informal assumptions of locality: the initial guess is not inverted, i.e. a(x1)> 0,and that the solver employed must follow a path which never jumps through thebarrier to inversion.We formalize these requirements in distortion optimization asminx{E(x) : a(x)≥ 0}. (3.15)Adding the constraint a(x) ≥ 0 now explicitly restricts our optimization to non-inverting deformations but otherwise leaves the desired solution unchanged. (SeeAppendix B for proof.)3.5.3 Iterating Away from CollapseWith problem statement Equation 3.15 in place, we now exploit it in curing thesearch direction from collapsing elements. At each iterate i, form the projectionminp{‖p+Di∇E(xi)‖22 : a(xi)+∇a(xi)T p≥ 0} (3.16)of the “uncured” predicted descent direction −Di∇E(xi) onto the subset satisfy-ing a linearization of the no-collapse condition. Satisfying Equation 3.16 exactlywould ensure that projected directions would not violate linear collapse conditionsand likewise preserve symmetry [118]. However, its exact solution is neither neces-60sary nor efficient. Instead, we construct an approximate solution to Equation 3.16as a filter that helps avoid collapse, preserves symmetry, and guarantees a low costfor computation for all descent steps.Strict convexity of the projection guarantees that a minimizer p∗ of Equa-tion 3.16 is given by the KKT5 conditions [15]p∗+Di∇E(xi)−∇a(xi)λ ∗ = 0, (3.17)0≤ λ ∗ ⊥ a(xi)+∇a(xi)T p∗ ≥ 0. (3.18)Setting Ci = ∇a(xi), Mi = ∇a(xi)T∇a(xi), and bi = a(xi), we form the Schur com-plement of the above system to arrive at the equivalent Linear ComplementarityProblem (LCP) [34]0≤ λ ∗ ⊥Miλ ∗+CTi pi+bi ≥ 0. (3.19)We then form a damped Jacobi splitting Mi =ω−1Si+(Mi−ω−1Si), with diagonalSi = diag(Mi) and damping parameter ω ∈ (0,1). This gives us an iterated LCPprocess ranging over iteration superscripts j,0≤ λ j+1 ⊥ ω−1Siλ j+1+Miλ j−ω−1Siλ j +CTi pi+bi ≥ 0. (3.20)3.5.4 Line Search FilteringEach iteration of the splitting Equation 3.20 simplifies to the damped projectedJacobi (DPJ) update6λ j+1← [λ j−ωS−1i (CTi (Ciλ j)+ ci)]+ , (3.21)with constant ci = CTi pi + bi. Here each of the m entries in λ j+1 can be updatedin parallel (unlike Gauss-Seidel iterations). As Mi is PSD this iteration processconverges to Equation 3.19 [34] and so to Equation 3.16. We do not seek a tight5Here and in the following λ = (λ1, ..,λm)T ∈Rm is the vector of Lagrange multipliers and x⊥ yis the complementarity condition xtyt = 0, ∀t.6We use the convention [·]+ = max[0, ·].61solution. Instead, we just want to be sure the worst blocks to line search are filteredaway. Thus we initialize with λ 1 = 0 to avoid unnecessary perturbation, use acoarse termination tolerance for DPJ (see below), and never apply more than amaximum of 20 DPJ iterations.At each DPJ iteration j we check for termination with an LCP specializedmeasure, the Fischer-Burmeister function [39] FB(λ j,Miλ j + ci) evaluated asFB(a,b) =√√√√ ∑k∈[1,m](ak +bk−√a2k +b2k)2. (3.22)As we initialize with λ 1 = 0, when pi is non-collapsing FB = 0, and thus no linesearch filtering iterations will be applied. Likewise, we stop iterations whenever theFB measure is roughly satisfied by either a relative error of < 10−3 or an absoluteerror < 10−6.Filtering thus applies a fixed maximum upper limit on computation and per-forms no iterations when not necessary. Upon termination of DPJ iterations, plug-ging our final λ into Equation 3.17 we obtain our update to form the line searchfiltered descent directionp`i = pi+Ciλ . (3.23)As illustrated in Figures 3.2 and 3.7 the filter’s small number of applied iterationscan make a dramatic difference in line search and so convergence.3.6 Termination CriteriaEvery iterative method for minimizing an objective function E(x) must incorpo-rate stopping criteria: when should an approximate solution be considered goodenough to stop and claim success? Clearly, in the usual case, where the actual localminimum value of E(x) is unknown, basing the test on the current value of E(xi) isfutile. As noted in Section 3.3.4, stopping when successive iterates are closer thansome tolerance is vulnerable to false positives (halting far from a solution), as isusing a fixed number of iterations. Although monitoring ‖∇E‖ is robust, each indi-vidual problem may need a different tolerance to define a satisfactory solution even62100 200 300 40034567(a)(b)0 10020030040034567SGDL-BFGSSL-BFGSAQPSGD and ConstraintsBCQN(a) AQP(b) SGD with filterwith filterIterationsFigure 3.7: Line search filtering. Bottom: We optimize a UV parameteriza-tion for a range of methods with the MIPS energy to consider line searchfiltering behavior, plotting energy (y-axis) against iteration counts. Justadding our barrier-aware line search filtering alone to SGD improves itsconvergence by well over an order of magnitude, and almost an orderof magnitude over AQP as well as plain L-BFGS and SL-BFGS. BCQNwith blending and line search filtering improves convergence even fur-ther. Top: a comparison between the embeddings and texture-mapsgained by AQP and barrier-aware filtered SGD at the 40th iterate.when normalized by number of vertices or total volume/area: see Figures 3.3 and3.11. We thus propose a new way to derive and construct an appropriate, roughlyproblem-independent, relative scale for a gradient-based measure for a stoppingcriterion.3.6.1 Characteristic Gradient NormAll energies we consider are summations of per-element energy densities W (·)computed from the deformation gradient Ft(x) and weights vt, in each element t, asper Equation 3.2. To simplify the following we can then evaluate energy densitieson the vectorized deformation gradient as W(vec(Ft))=W (Gtx), where Gt is the63linear gradient operator for element t. The full energy gradient is then∇E(x) = ∑t∈TvtGTt ∇W(Gtx). (3.24)We wish to generate a characteristic value we can compare this gradient to mean-ingfully, with the same dimensions; we will do this with each component of theabove summation separately.First observe that the deformation gradient, Ft , the argument to W , is dimen-sionless and therefore ∇W has the same dimensions as W , and even as the elementHessian ∇2W . For the simplest quadratic energy densities, this Hessian has theattractive property of being constant; we thus choose to use the 2-norm of theHessian, evaluated about the deformation gradient at rest (Ft = I), to get a repre-sentative value for ∇W :〈W 〉= ‖∇2W (I)‖2. (3.25)Second, note that the ith part of Gt for a triangle (respectively tetrahedra) t con-taining vertex i will attain its maximum value for fields which are constant alongthe opposing edge (triangle) and that value will be the reciprocal of the altitude. Upto a factor of 2 (3), this is the length (area) of the opposing edge (triangle) dividedby the rest area (volume), of the element, i.e. vt. Summing over all incident ele-ments, weighted by vt, we arrive at a characteristic value for vertex i of `i equallingthe perimeter (surface area) of the one-ring of vertex i. We compute this value `i,for all vertices, giving us the vector `(T ) = (`1, ..., `n)T ∈Rn, with one scalar entryper vertex.The product of our energy and mesh values together form the characteristicvalue for the norm of the gradient〈W 〉‖`(T )‖, (3.26)where we take the same vector norm as that with which we evaluate ‖∇E(x)‖; weuse the 2-norm in all our experiments. For all methods we stop iterating when‖∇E(x)‖ ≤ ε〈W 〉‖`(T )‖, (3.27)64given a dimensionless tolerance ε from the user, which is now essentially mesh-and energy-independent. See Figures 3.3, 3.10 and 3.11 as well as our experimentalanalysis in Section 3.8 for evaluation.3.7 The BCQN AlgorthimAlgorithm 2 contains our full BCQN algorithm in pseudocode. The dominant cost,for both memory and runtime, is the Laplacian solve embedded in the applicationof D, which again is not stored as a single matrix, but rather is a linear transforma-tion involving a few sparse triangular solves with the Laplacian’s Cholesky factorand outer-product updates with a small fixed number of L-BFGS history vectors.Recall that we separately solve for each coordinate with a scalar Laplacian, not us-ing a larger vector Laplacian on all coordinates simultaneously; this also exposessome trivial parallelism. Apart from the Laplacian, all steps are either linear (dot-products, vector updates, gradient evaluations, etc.) or typically sublinear (DPJassembly and iterations, which only operate on the small number of collapsingtriangles, and again are easily parallelized).As Lipton et al. proved [77], the lower bounds for Cholesky factorization ona two-dimensional mesh problem with n degrees of freedom are O(n logn) spaceand O(n3/2) sequential time, and in three-dimensional problems where vertex sep-arators are at least O(n2/3), their Theorem 10 shows the lower bounds are O(n4/3)space and O(n2) sequential time. On moderate size problems running on currentcomputers, the cost to transfer memory tends to dominate arithmetic, so the spacebound is more critical until very large problem sizes are reached.3.7.1 Comparison with Other AlgorithmsThe per-iterate performance profile of AQP is most similar to BCQN: it too is dom-inated by a Laplacian solve. The only difference is the extra linear and sublinearwork which BCQN does for the quasi-Newton update and the barrier-aware fil-tering; even on small problems, this overhead is usually well under half the timeBCQN spends, and as the next section will show, the improved convergence prop-erties of BCQN render it faster.The second-order methods we compare against, PN and CM, as well as the65Algorithm 2 Blended Cured Quasi-Newton (BCQN)Given: x1, E, ε , TInitialize and Precompute:r = ε〈W 〉‖`(T )‖ // Characteristic termination value (§3.6)L, D← L−1 // Initialize blend model (§3.4)g1 = ∇E(x1), i = 1while ‖gi‖> r do // Termination criteria (§3.6)p←−Dgi // Precondition gradient (§3.4)// Assemble for DPJ iterations (§3.5):C← ∇a(xi)M←CTC, c←CT p+a(xi)S← diag(M)−1, λ ← 0fb← FB(λ ,Mλ + c) // LCP residual (Equation (3.22) in §3.5)for j = 1 to 20 // Line-search filtering (§3.5)if fb< 10−6 then break end iffb← fbnextλ ← [λ − 12 S(CT (Cλ )+ c)]+ // Parallel project (§3.5)fbnext← FB(λ ,Mλ + c)if |fb− fbnext|/fb< 10−3 then break end ifend forp`← p+Cλ // Line-search filtered search direction (§3.5)α ← LineSearch(xi, p`,E) // Line search (§3.4)xi+1 = xi+α p` // Descent step (§3.4)gi+1 = ∇E(xi+1)D← Blend(D,L,xi+1,xi,gi+1,gi) // BCQN blending update (§3.4)i← i+1end whilemore approximate proxy methods, SLIM and AKAP, all use a fuller stencil whichcouples coordinates. The same asymptotics for Cholesky apply, but whereas AQPand BCQN can solve a scalar n× n Laplacian d times (once for each coordinate,independently), these other methods must solve a single denser nd× nd matrix,with d2 times more nonzeros: see Figure 3.8. Moreover, for all these methods theproxy matrix changes at each iteration and must be refactored, adding substantiallyto the cost: factorizations are significantly slower than backsolves.66LaplacianHessianFill (AMD Order)XYZX Y ZXYZX Y ZFigure 3.8: Sparsity Differences in Proxies. Left: The scalar Laplacian(top) is smaller and sparser than the Hessian and its approximations(bottom) used in CM, PN, SLIM and AKAP. Right: This results in amuch cheaper factorization and solve for the Laplacian; it is appliedin both BCQN and AQP independently to each coordinate where onlya one-time factorization precompute is required; CM, PN, SLIM andAKAP require factorization at each iterate.3.8 Evaluation3.8.1 ImplementationWe implemented a common test-harness code to enable the consistent evaluationof the comparative performance and convergence behavior of SGD, PN, CM, AQP,L-BFGS and BCQN across a range of energies and distortion optimization tasksincluding parameterization as well as 2D and 3D deformations, where these meth-ods allow. For AQP this extends the number of energies it can be tested with, while67more generally providing a consistent environment for evaluating all methods. Wehope that this code will also help support the future evaluation and development ofnew methods for distortion optimization.The main body of the test code is in MATLAB to support rapid prototyping. Alllinear system solves are performed with MATLAB’s native calls to CHOLMOD[29] with additional computational-heavy modules, primarily common energy, gra-dient and iterative LCP evaluations, implemented in C++. As linear solves are thebottleneck in all methods covered here, an additional speed-up to all methods ispossible with Pardiso [101, 102] in place of CHOLMOD; however, as discussed inSection 3.8.4 this does not change the relative merits of the methods, and wouldadd an additional external dependency to the test code. For verification we alsoconfirm that iterations in the test-harness AQP and CM implementations match theofficial AQP [70] and CM [115] codes.All experiments were timed on a four-core Intel 3.50GHz CPU. We have par-allelized the damped Jacobi LCP iterations with Intel TBB; with more cores theoverhead reported below for LCP iterations is expected to diminish rapidly. For allUV parameterization problems we follow Kovalsky et al. [70] and compute an ini-tial, locally injective embedding via a single linear solve with the cotan Laplacian;if it fails (is not locally injective), we then fall back to plain Tutte, so that robustnessis maintained. For all constrained deformation examples, with the exception of theArmadillo in Figure 3.14, we begin with an initially injective mapping. For theArmadillo deformation example only, we apply the LBD method [69] to create arough, locally injective initialization from the initial constrained non-injective de-formation and pass this to both compared methods. To enforce Dirichlet boundaryconditions, i.e. positional constraints, we use a standard subspace projection [93],i.e. removing those degrees of freedom from the problem. For constrained prob-lems we apply the standard approach of projecting gradients to the null-space ofconstraints: a stationary point is then reached only if the (projected) gradient van-ishes, just as in the unconstrained setting. When line search is employed we firstfind a maximal non-inverting step size with Smith and Schaefer’s method [120],followed by standard line search with Armijo and curvature conditions.68krE(xi)kCharacteristickE(xi)E(xi1)k/kE(xi)k1e-1 1e-2 1e-3 1e-4 1e-5 1e-6?? 9e-61e-1 1e-3 8e-9 1e-9 3e-10Swirl OptimizationFigure 3.9: Termination criteria comparison. Left to right: We find keypoints in the sequential progress of the optimized mesh in the Swirloptimization (ISO energy) example at regular intervals of 10× decreasein our characteristic norm. We compare with the relative error measuresat these same points.3.8.2 TerminationTo evaluate termination criteria behavior we first instrumented two geometry op-timization stress-test examples: the Swirl deformation [28] and the Hilbert curveUV parametrization [120]. We run both examples to convergence (ε = 10−6 usingour characteristic gradient) reaching the final target shapes for each. Within theseoptimizations we record the 2-norm of gradient, the vertex-normalized 2-norm ofgradient, the relative error measure [70, 115] and our characteristic gradient normfor all iterations.Figure 3.9 shows the Swirl mesh obtained during BCQN iteration at regularintervals of 10× decrease in our characteristic norm. Observe that they correspondto natural points of progress. For comparison we also provide the correspondingrelative error measures, which varies much less steadily.In Figure 3.10 we compare termination criteria more closely for a parametriza-tion problem, the Hilbert curve example. We plot our characteristic gradient norm(blue) and the relative energy error [70, 115] (orange) as BCQN proceeds. Notethat the characteristic gradient norm provides consistent decrease corresponding toimproved shapes and so provides a practical measure of improvement. The localerror in energy, on the other hand, varies greatly, making it impossible to judgehow much global progress has been made towards the optimum.Figure 3.11 illustrates consistency across changing tolerance values, mesh res-olutions, and scales. example. We show the iterates at measures 10−3, 10−4and691 5 10 50 100 500 100010-60.00111000106OursYaronFigure 3.10: Measuring improvement. Solving a UV parametrization of theHilbert curve with BCQN, we plot our characteristic gradient norm inblue and the relative energy error in orange as the method proceeds,on a logarithmic scale. Iterates are shown at decreases in the char-acteristic gradient norm by factors of 10, illustrating its efficacy as aglobal measure of progress, while the relative energy error measuresonly local changes with little overall trend.10−5 for both our characteristic gradient norm and the raw gradient norm, formeshes with varying refinement and varying dimension (rescaling coordinates bya large factor). Similar to Figure 3.3 comparing the vertex-normalized gradientnorm, there are large disparities for the raw gradient norm, but our characteristicgradient norm is consistent.Tolerances The Swirl and Hilbert curve examples are both extreme stress teststhat require passing through low curvature regions to transition from unfolding tofolding; see e.g., Figure 3.9 above and our videos. For these extreme tests we useda tolerance of 10−6 for our characteristic gradient norm to consistently reach thefinal target shape. However, for most practical distortion optimization tasks sucha tolerance is excessively precise. In experiments across a wide range of energies70CoarseFineScaled1e-3 1e-4 1e-51e-3 1e-4 1e-5CharacteristicTolerancekrE(xi)k krE(xi)kFigure 3.11: Termination criteria comparison across mesh refinement andscale. Left and right: we show the Swirl optimization when ourcharacteristic norm (left) and the standard gradient norm (right) reach10−3,10−4 and 10−5. Top to bottom: the rows show optimizationwith a coarse mesh, a fine mesh, and the same fine mesh uniformlyscaled in dimension by 100×. Note the consistency across mesh reso-lution and scaling for our characteristic norm and the disparity acrossthe standard gradient norm.and UV parametrization, 2D and 3D deformation tasks, including those detailedbelow, we found that ‖∇E(x)‖ ≤ 10−3〈W 〉‖`(V,T )‖ consistently obtained good-looking solutions with essentially no further visible (or energy value) improvementpossible. We argue this is a sensible default except in pathological examples. Forall examples discussed here and below, with the exception of the Swirl and theHilbert curve tests, we thus use ε = 10−3 for testing termination.3.8.3 Newton-type MethodsWhile Newton’s method, on its own, handles energies like ARAP reasonably well[25] it is insufficient for strongly nonconvex energies: modification of the Hessianis required [93, 115]. Here we examine the convergence, performance and scal-ability of PN [134], a general-purpose modification for nonconvex energies, andCM [115], a more recent convex majorizer currently restricted to 2D problems and71Vertices TrianglesBCQN AQP Projected Newton Composite Majorization1.9K 3.1K 95 0.98 0.02M 126 1.25 0.02M 19 1.01 0.07M 19 0.84 0.07M3.5K 6.3K 86 1.44 0.05M 148 2.30 0.05M 17 1.86 0.19M 19 1.57 0.19M6.6K 12.5K 78 2.18 0.11M 153 4.67 0.11M 27 6.01 0.45M 19 3.26 0.45M12.9K 25.0K 118 6.51 0.27M 197 10.27 0.27M 26 12.81 1.02M 19 7.82 1.02M25.4K 50.0K 110 11.21 0.60M 223 25.23 0.60M 19 19.83 2.38M 20 18.02 2.38M50.4K 100.0K 101 22.41 1.41M 177 39.41 1.41M 25 56.10 5.56M 20 40.89 5.56M100.4K 200.0K 104 48.10 3.24M 249 115.06 3.24M 19 95.58 12.55M 20 99.43 12.55M197.9K 394.6K 112 111.08 6.14M 225 213.04 6.14M 28 290.26 24.13M 20 208.83 24.13M435.5K 869.2K 119 355.44 15.37M 261 758.92 15.37M 29 757.01 56.94M 20 593.88 56.94M880.3K 1,758.1K 138 758.12 36.59M 341 1,750.50 36.59M 36 2,160.40 125.89M 21 1,306.64 125.89M1,650.4K 3,297.5K 151 1,532.10 75.05M 377 3,767.20 75.05M 27 4,461.30 395.71M 21 3,500.21 395.71M3,221.7K 6,438.7K 107 3,584.80 176.55M 354 7,304.90 176.55M 44 40,788.00 709.36M 21 13,024.00 709.36M6,386.2K 12,765.6K 145 19,209.00 347.87M 317 42,895.00 347.87M * 1,385.60M * 1,385.60M11,969.0K 23,928.4K 138 54,550.00 727.61M 452 169,460.00 727.61M * 3,771.30M * 3,771.30MIteration Timing(s) Fill-in Iteration Timing(s) Fill-in Iteration Timing(s) Fill-in Iteration Timing(s) Fill-inFigure 3.12: UV Parameterization Scaling, Timing and Sparsity. Perfor-mance statistics and memory use for increasing mesh sizes up to 23.9Mtriangles, comparing BCQN with AQP, PN and CM. For the GorillaUV parametrization with ISO energy we repeatedly double the meshresolution and, for each method, report number of iterations to con-vergence (characteristic norm < 10−3), wall-clock time (seconds) toconvergence, and the nonzero fill-in for the linear systems solved byeach method. We use * to indicate out-of-memory failure for matrixfactorization; see §3.8 for discussion. Also note that stencils for CMand PN are identical (differing only by actual entries) while AQP andBCQN both solve with the same smaller scalar Laplacian.a trio of energies (ISO, Symmetric ARAP and NH), and compare them with AQPand BCQN. For the 2D parameterization problems in Figure 3.12 we can compareall four methods while for the 3D deformation problems in Figures 3.13 and 3.14CM is not applicable.Figures 3.12, 3.13, and 3.14 examine the scaling behavior of the various meth-ods under mesh refinement, for 2D parameterization and 3D deformation. TheNewton-type methods PN and CM (when applicable) maintain low iteration countsthat only grow slowly with increasing mesh size; from the outset BCQN and AQPrequire more iterations, though the iteration count also grows slowly for BCQN.72Nonetheless, BCQN is the fastest across all scales in each test as its overall costper iteration remains much lower. BCQN iterations require no re-factorizations(which scales poorly, particularly in 3D, as discussed in Section 3.7) and onlysolves smaller and sparser scalar Laplacian problems per coordinate compared tothe larger and denser system of CM and PN. This advantage for BCQN only in-creases as problem size grows; indeed, for the largest problems BCQN succeededwhere CM and PN ran out of memory for factorization.3.8.4 A Note on Solving Proxies and PardisoRecent methods including CM have taken advantage of the efficiencies and opti-mizations provided by the Pardiso solver. While this can improve runtime of thefactorization and backsolves by a constant factor, it cannot change the asymptoticlower bounds on complexity; the sparse matrix orderings in both CHOLMOD andPardiso already appear to achieve the bound on typical mesh problems. In tests onour computer, across a large range of scales in two and three dimensions, we foundPardiso was occasionally slower than CHOLMOD but usually 1.4 to 3 times faster,and at most 8.1 times faster (for backsolving with a 3D scalar Laplacian).Individual iterates of AQP have the same overall efficiency as BCQN (domi-nated by the linear solves); switching to Pardiso leaves the relative performanceof the two methods unchanged. While CM and PN are even more dependent onthe efficiency of the linear solver, due to more costly refactorization each step, thesame speed-ups possible with Pardiso also apply to BCQN, so again there is nosignificant change in relative performance between the methods.3.8.5 First-order MethodsAmong existing first-order methods for distortion optimization AQP has so farshown best efficiency [70] with improved convergence over SGD as well as stan-dard L-BFGS. Likewise, as we see in Figures 3.12, 3.13 and 3.14, when we scaleto increasingly larger problems AQP will dominate over Newton-type methods andso potentially offers the promise of reliability across applications. Finally, BQCNperforms a small, fixed amount of extra work per-iteration in the line-search filterand quasi-Newton update. Thus in Figures 3.12, 3.13, 3.15, 3.16 and 3.17 we com-73Vertices Tetrahedra BCQN AQP Projected Newton1.8K 6.3K 62 4.22 262 15.97 41 22.01 0.61M 0.03M3.0K 10.8K 63 7.02 311 30.56 42 39.61 1.35M 0.07M6.6K 25.6K 76 18.27 342 79.51 43 97.21 4.88M 0.27M16.0K 66.6K 130 74.90 408 226.02 44 282.64 21.04M 1.09M31.7K 137.2K 151 168.59 475 517.32 44 644.97 65.06M 3.29M65.3K 291.6K 178 429.48 546 1,346.30 44 1,673.90 214.44M 10.44M133.1K 608.4K 274 1,391.70 538 2,640.00 44 15,650.00 695.33M 31.25M261.9K 1,219.5K 291 2,805.10 685 7,300.00 45 76,711.00 1,135.20M 87.26M535.7K 2,532.7K 352 7,930.60 825 28,500.00 45 1,104,480.00 3,078.80M 270.73M1,633.5K 7,873.2K 392 20,908.00 919 88,215.00 * 9,672.17M 689.42MIteration Timing(s) Iteration Timing(s) Iteration Timing(s) Fill-inLaplacian Fill-inFigure 3.13: Three-Dimensional Deformation Scaling, Timing and Spar-sity. Performance statistics and memory use for increasing mesh sizesup to 7.8M tetrahedra, comparing BCQN with AQP and PN. (CM doesnot extend to 3D.) We initialize a bar with a straight rest shape to startin a tightly twisted shape, constraining both ends to stay fixed and thenoptimize over increasing resolutions. For each method we report num-ber of iterations to convergence (characteristic norm< 10−3) and wall-clock time (seconds) to convergence, as well as the nonzero fill-in forthe linear system solved by CM as compared with the Laplacian em-ployed by both BCQN and AQP. We use * to indicate out of memoryfor the computation on our test system; see §3.8 for discussion.pare AQP and BCQN over a range of practical distortion optimization applications,UV-parameterization, 2D deformation, and 3D deformation, with nonconvex ener-gies from geometry processing and physics. Throughout we note three key featuresdistinguishing BCQN as follows.Reliability and robustness: AQP will fail to converge in some cases, see e.g.Figure 3.4, while BCQN reliably converges. In our testing AQP fails to convergein over 40% of our tests with nonconvex energies; see e.g. Figures 3.15 and 3.16.This behavior is duplicated in our test-harness code and AQP’s reference imple-mentation.Convergence speed: when AQP is able to converge, BCQN consistently pro-74Locally Injective Initializer Local MinimizerRest ShapeApplied ConstraintsFigure 3.14: Armadillo Deformation Test. We compare three-dimensionaldeformation optimizations of a 1.5M element tetrahedral mesh of theT-pose armadillo with BCQN and PN. We constrain the armadillo’sfeet to rest position, its hands to touch the ground and use the LBDmethod [69] to create a locally injective initialization for the solversfrom the starting constrained configuration. Here BCQN requires 393iterations to converge while PN converges in just 9. However, asBCQN remains cheaper and more scalable per iterate it takes only4,148 seconds, while PN spends 13,447 seconds in our test system.See Figure 3.8 for the comparison of fill-in and size.vides faster convergence rates for nonconvex energies. In our experiments conver-gence rates range up to well over 10X with respect to AQP.Performance: BCQN is efficient. When AQP is able to converge, BCQN re-mains fast with up to a well over 7X speedup over AQP on nonconvex energies.75Example Vertices Triangles EnergyBCQN AQPIteration Timing Iteration TimingBull 17.9K 34.5KISO 169 11.82 † †MIPS 200 24.08 † †CONF 425 50.06 1,085 74.45Camel 40.2K 78.1KISO 412 92.60 593 114.89MIPS 3,344 821.57 66,479 12,281.00CONF 9,196 2,067.60 81,100 15,712.00Dino 24.6K 47.9KISO 162 30.43 950 98.66MIPS 283 54.13 † †CONF 156 29.53 883 98.02Isis 188.1K 374.3KISO 347 404.63 353 366.01MIPS 204 298.55 † †CONF 104 166.65 611 612.98Cow 3.1K 5.8KISO 78 1.14 † †MIPS 239 3.29 323 3.98CONF 59 1.05 108 1.23Horse 20.6K 39.6KISO 104 12.93 643 57.40MIPS 387 47.73 † †CONF 519 68.22 5,552 523.97Figure 3.15: UV parameterization. Top row: 3D meshes for UVparametrization with ISO, MIPS, and CONF distortion energies. Mid-dle two rows: converged maps and texturing from BCQN on ISO ex-amples. Bottom: for each method / problem pair we report number ofiterations to convergence (characteristic norm < 10−3) and wall-clocktime (seconds) to convergence. We use † to indicate when AQP doesnot converge; see §3.8.5.763.8.6 Across-the-board ComparisonsHere we compare the performance and memory usage of BCQN with best-in-classdistortion optimization methods across the board: AQP, PN and CM for both 2Dparameterization and 3D deformation tasks. Results are summarized in Figures3.12, 3.13 and 3.14. Note that CM does not extend to 3D.In Figures 3.12 and 3.13 we examine the scaling of AQP, PN, CM and BCQNto larger meshes and thus to larger problem sizes in both 2D parametrization (upto 23.9M triangles) and 3D deformation (up to 7.8M tetrahedra). As noted above:from the outset, BCQN requires more iterations than CM and PN; however, BCQN’soverall low cost per iteration makes it faster in performance across problem sizeswhen compared to both CM and PN. We then note that AQP, on the other hand,has slower convergence and so, at smaller sizes it often does not compete with CMand PN. However, once we reach larger mesh problems, e.g. ∼≥ 6M triangles inFigure 3.12, the cost of factorization and backsolve of the denser linear systems ofCM and PN becomes significant so that even AQP’s slower convergence results inimprovement. This is the intended domain for which first-order methods are de-signed but here too, as we see in Figure 3.12, BCQN continues to outperform bothAQP as well as CM and PN across all scales.3.9 ConclusionIn this work we have taken new steps to both advance the state of the art for op-timizing challenging nonconvex deformation energies and to better evaluate newand improved methods as they are subsequently developed. Looking forward theseminimization tasks are likely to remain fundamental bottlenecks in practical codesand so advancement here is critical. Our three primary contributions together formthe BCQN algorithm which pushes current limits in optimizing distortion forwardin terms of speed, reliability, and automatability. At the same time, looking aheadwe also expect that each contribution individually should lead to even more im-provements in the near future.77Example Vertices Triangles EnergyBCQN AQPIteration Timing Iteration TimingSwirl 0.5K 0.8KISO 259 0.32 † †MIPS 55 0.13 220 0.22NH 202 0.43 † †Shear 5K 9.7KISO 192 4.56 † †MIPS 435 9.25 † †NH 710 15.07 † †Gecko 0.8K 1.2KISO 65 0.17 † †MIPS 87 0.22 † †NH 62 0.16 † †Bend 0.5K 0.8KISO 25 0.05 † †MIPS 23 0.05 † †NH 18 0.04 † †Figure 3.16: Two-Dimensional Deformation. Top: initial conditions andvertex constraints (blue points) for deformation problems minimizingISO, MIPS, and NH deformation energies. Middle: converged so-lutions from BCQN on ISO examples. Bottom: for each method /problem pair we report number of iterations to convergence (charac-teristic norm < 10−3) and wall-clock time (seconds) to convergence.We use † to indicate when AQP does not converge; see §3.8.5.78Example Vertices Tetrahedra EnergyBCQN AQPIteration Timing Iteration TimingDancer 15.6K 54.9KISO 38 23.49 55 33.56MR 56 43.32 87 68.49Botijo 15.9K 57.0KISO 11 8.54 49 45.20MR 15 12.35 65 59.95Dino 16.1K 58.1KISO 15 11.77 17 9.31MR 15 13.84 40 23.37Homer 16.5K 61.2KISO 18 13.94 38 23.56MR 11 8.66 72 45.64Figure 3.17: Three-Dimensional Deformation. Top: initial conditions forvertex-constrained deformation problems minimizing ISO and MRdeformation energies. Middle: converged solutions satisfying con-straints from BCQN on MR examples. Bottom: for each method /problem pair we report number of iterations to convergence (charac-teristic norm < 10−3) and wall-clock time (seconds) to convergence.793.9.1 Limitations and Future WorkWhile our focus is on recent challenging nonconvex energies not addressed bythe popular local-global framework, similar to AQP we have observed significantspeedup for energies that are convex in singular values, e.g. ARAP, as well. Cur-rently, in comparing AQP and BCQN on the same set of 2D and 3D tasks withthe ARAP energy we observe a generally modest improvement in convergence, upto a little over 4×, which is somewhat balanced by the small additional overheadof BCQN iterations. Note for energies like ARAP there is no barrier, hence noneed for our line search filtering, but other opportunities for improvement may beavailable in future research.While our current blending model works well in our extensive testing, it is em-pirically constructed; it is in no sense proven optimal. We believe that further anal-ysis, better understanding and additional improvements in quasi-Newton blendingare all exciting and promising avenues of future investigation.Our focus has been on unconstrained optimization, such as problems where lin-ear constraints (e.g. pinned vertices) have already been removed from the degreesof freedom. As stated earlier, unconstrained optimization is also often a crucialinner step in constrained optimization methods such as Augmented Lagrangian orInterior Point. However, we caution that while our filtering step, based on lin-earized constraints to prevent element inversion, is a sound and effective means toimprove the subsequent line search, it is not a general-purpose tool to handle evenjust linear constraints: it relies crucially on the fact that the no-inversion constraintswe add should already be satisfied by the unconstrained solution, i.e. ultimately areinactive. We leave the question of incorporating active constraints directly into fil-tering and blending BFGS for future work and, as a starting point, observe thatlinear constraints can be directly added with a KKT formulation.We also note that while we have focused here on optimizing distortion energiesdefined on meshes, there is a wide range of critical optimization problems that takea similar general structure: minimizing separable nonlinear energies on graphs.Further extensions are thus exciting directions of ongoing investigation.Finally, one limitation of our solver is that its performance scales with thevolumetric mesh resolution. To account for this issue, we propose a boundary inte-80gral technique in the next chapter formulating the optimization problem underlyingquasistatic elasticity over the boundary mesh, eliminating the need for a volumetricmesh and with it the related performance scaling trouble.81Chapter 4Simulating Rigid Body Fracturewith Surface MeshesHigh-quality rigid body fracture has been a popular research topic in computergraphics, with the purpose of increasing realism and plausibility in computer-generated animation. Fracturing effects such as shattering glass, destruction ofa cement wall or blowing up a space battleship greatly affect the immersive userexperience in computer games, virtual reality and film industry. Roughly speaking,current state-of-the-art approaches can be categorized into two groups: (1) crackpropagation and (2) volume decomposition. The first category uses a volume dis-cretization to solve for the underlying elastic dynamics, and applies some crackgrowth criteria, such as principal stress hypothesis (Rankine hypothesis), criticalstrain energy release rate (Griffith criterion) and minimum strain energy density(Beltrami hypothesis), to propagate fractures inside the material. These methodsare generally physically realistic, but entail significant computational cost due toeither the poor scaling behavior of the physical system or the complexity of the un-structured volumetric grid’s remeshing operation. The second category starts withspatial decomposition of the object instantly, instead of waiting for the crack prop-agation process to create new fractured pieces. Such a simplified strategy gives riseto computational efficiency and flexible control of fragment generation, which isessential and attractive to graphics experts, practitioners and users. However, thepotentially unpredictable physically inconsistent behavior caused by this simplifi-82cation might require significant post-processing.Our goal is to design a new rigid body fracture algorithm that is physicallyplausible and computationally efficient. We present a surface-only method to han-dle both physical computation and fracture evolution by formulating the quasi-static elasticity problem as an indirect boundary integral equation and using ex-plicit mesh-based interface tracking technique (see Figure 4.1). By making use ofthe indirect formulation, we eliminate the non-trivial null space in the system offree-flying objects, which is an inherent problem for the direct formulation. Elas-tic stress tensors, which are the crucial ingredient for our fracture initiation andfracture propagation method, can be thus written in a consistent analytical form.We initiate and propagate crack interface in the spirit of the Rankine hypothesis.Each crack front will be generated and evolved under the principal stress flow inthe object’s material space explicitly and independently. We use existing robustmesh-based surface tracking methods, but incorporate significant modifications inorder to adapt the tracking algorithm to our fracture evolution problem. The ex-plicit approach not only avoids the need for a signed distance volume grid but alsopreserves tiny surface geometrical features. Another benefit is that our method canbe easily combined with any popular rigid body dynamics engine by using an ex-plicit surface representation. In our implementation, we use a velocity level linearcomplementary programming (LCP) rigid body solver. The physical equation us-ing our indirect formulation is well-conditioned, and discretization of the integralequation leads to a dense system, which is small relative to the the sparse systemthat arises in the volumetric approach, but it nonetheless requires a sophisticatedsolution approach. To that end, we modify a kernel-independent Fast MultipoleMethod and combine it with a BiCG-STAB iterative solver, to obtain a linear run-time in the number of surface triangles.4.1 IntroductionIn summary, our core technical contributions are as follows:• we have developed the first rigid body fracturing physical simulation thatonly requires surface representation and surface mesh operations;• we have designed and implemented a new well-conditioned indirect bound-83Figure 4.1: A glass goblet broken into multiple shards (left), with corre-sponding fracture propagation shown in the material space (right). Onlythe surface mesh is used for both the physical computation and fracturepropagation.84ary integral formulation for the linear quasi-static elasticity problem;• we offer the first explicit mesh-based surface tracking algorithm for the rigidbody brittle fracture problem.We also highlight the flexibility of the kernel-independent Fast Multipole Methoddesigned for the linear time n-body fast summation, which to the best of our knowl-edge has not been previously used in computer graphics models and applications.4.2 Related WorkFracture Simulation. Since the pioneering work of Terzopoulos and Fleischer[135] and Norton et al. [94], fracture modeling has been an important part of phys-ically based animation and interactive virtual environment. In their seminal work,O’Brien and his collaborators [95, 96] simulated brittle and ductile fracture bymaking use of the finite element method (FEM). In order to deal with the elastody-namic problem, researchers used a corotational formulation to handle large defor-mation and fracturing effects [55, 88, 99]. One of the most challenging problems inFEM-based fracture simulation arises from the topology change of the volumetricmesh. Mesh elements cut by the crack front must be dynamically remeshed so as tohave the simulation domain conform with the fracture surface [95, 143]. A naı¨vetreatment, like deleting elements or aligning the crack plane with the elements’boundary, may lead to a jagged result. Molino et al. [85] and Sifakis et al. [116]introduced a virtual node algorithm to decouple material domain and simulationdomain by duplicating elements necessary for modeling the topological change.To avoid volumetric mesh updates, Pauly et al. [100] and Steinemann et al. [125]developed meshless fracturing methods as a composition of mesh free simulationand crack front tracking. Unlike these fracture propagation approaches based onsurface reconstruction or crack synthesis, we develop a new crack front trackingalgorithm using simple triangle mesh operations similar to explicit mesh-basedsurface tracking techniques [21, 35]. Moreover, by formulating the linear quasi-static stress analysis [7, 89, 154] in our indirect boundary integral formulation, weplace all degrees of freedom on the surface only and get rid of the visibility graphas required by previous meshless methods. Other interesting fracturing work in-85cludes thin sheet problem [24, 61, 103], Voronoi diagram methods [7, 90, 126]and material strength field synthesis [30].Boundary Integral Equations. Boundary integral equations (BIE) are refor-mulations of boundary value problems for partial differential equations (PDEs).They have been present in the mathematics literature for almost two centuries, andengineers have made much use of them for solving a wide range of real-worldproblems. Their promise, motivating the present work, is to significantly reducethe number of degrees of freedom compared to volumetric methods, and simplifythe geometric part of the problem to surface meshes only. For example, to solveelasticity in a region with an n-element surface mesh, a BIE solver acceleratedby FMM or similar strategies takes O(n) time, and needs nothing more than thesurface mesh. By comparison, a uniform tetrahedral mesh conforming to the sur-face would contain O(n1.5) elements, and the common solver choice of a sparseCholesky factorization would take O(n2) space and O(n3) time.The use of BIE methods in physical animation dates back to the work of Jamesand Pai [56], who formulated the elastostatic problem as direct boundary integral.Such a formulation is suitable for elliptic PDEs with Dirichlet boundary condi-tions, but it has a nontrivial null space with pure Neumann boundary condition,for example a free flying rigid object. Researchers addressed this problem withinthe FEM framework and overcame it by anchoring degrees of freedom or project-ing out the null space algebraically [7, 89, 154]. Our own approach is to tacklethis ill-conditioned PDE by exploring an indirect boundary integral formulationwhich has previously been applied to vortex sheets [23] and ocean waves [64].A primary characteristic of our method is that it involves solving a second kindFredholm integral equation, which is guaranteed to be well-conditioned. More de-tails are provided in Section 4.4.2. Other graphics applications of BIE include thevortex method [43, 98, 141], time harmonic sound rendering [58] and diffusioncurve imaging [127]. BIE has also been widely explored in mechanic field, such aselstostatic problem, elastodynamics and fracture analysis, using various discretiza-tion approaches including collocation scheme [16], Bubnov Galerkin [40], PetrovGalerkin [117], etc. Approaches considering elastic wave include time-domainboundary element formulation [84], dual reciprocity method [105] as well as time-harmonic solution [109]. Since dynamic approach involves special treatment on86crack tip during fracture propagation to keep numerical accuracy [5], we insteadlook at the problem from a quasistatic point of view. We also refer interested read-ers to [31] for more discussions on singular integral quadrature.Crack Front Tracking. Conventional FEM-based Lagrangian fracture meth-ods use volumetric mesh cutting algorithms to keep track of fracture propagation.See [148] for a comprehensive and detailed review of previous cutting algorithmsin the literature. Our fracture propagation algorithm is mostly related to the mesh-based surface tracking techniques [21, 35] which have been successfully appliedto multiphase flow. We make several essential modifications to the mesh opera-tion strategy of the LosTopos library [35] in order to resolve the particular physicalbehaviour of crack initiation and propagation. Our mesh operation strategy doesdepend on the crack opening criterion. We use the Rankine condition, because itseems to be much more effective than the Griffith energy type approach [50] in pre-serving volume. The Rankine approach provides us with the flexibility to controlcrack open thickness, and it yields a physically convincing simulation result.Fast Integral Equation Solver. Compared with FEM or the finite differencemethod (FDM) for PDEs, the boundary integral formulation has the advantage ofhaving fewer degrees of freedom. However, the fact that unknown variables in BIEare all coupled leads to dense system after discretization, which hinders the di-rect application of a powerful iterative linear system solver. By using the low-rankproperty of n-body interaction with non-oscillatory kernel, algorithms of O(n) orO(nlog n) complexity for matrix-vector multiplication have been developed since1980’s, including particle-particle/particle-mesh (P3M) [52], tree code [9], panelclustering [46] and fast multipole method (FMM) [44]. A second class of meth-ods are wavelet compression techniques [4], which lead to sparse and asymptot-ically well-conditioned approximations of the coefficient matrix. Applications ofthe fast iterative dense system solver in computer graphics include deformable ma-terial [57], radiosity rendering [47], fluid simulation [23, 64] and sound rendering[154]. Instead of applying a spherical multipole expansion to the kernel functionwhich is a second order tensor in the elastostatic case, we introduce a radial basisfunction (RBF) style based kernel independent FMM in Section 4.5 by simplifyingthe implementation of Ying et al. [153]. By utilizing the RBF field approximation,we are able to write the multipole expansion, multipole to multipole (M2M), multi-87pole to local (M2L), local to local (L2L) and local expansion in a simple consistentform.Algorithm 3 Rigid Body Fracture Simulation Loop1: solve rigid body dynamics [Figure 4.2a]2: estimate rigid body contact force [Figure 4.2b]3: solve layer potential for elastostatic problem [Figure 4.2c]4: apply stress analysis on surface nodes [Figure 4.2d]5: add local maximum nodes to fracture node listN6: set fracture line listL ←Crack Initiation(N )7: for each line l inL8: label all nodes, except the ends, of l as active9: end for10: whileL 6= /0 do11: apply stress analysis on active crack nodes ofL [Figure 4.2d]12: label nodes as inactive if they fail on Rankine condition13: assign propagation velocity to active crack nodes [Equation 4.9]14: Crack Propagation(L , sub timestep) [Figure 4.2e]15: end while16: update rigid body list if new fragments are created4.3 Algorithm OverviewOur rigid body fracture simulation method is closely related to the work by Mu¨lleret al. [89], Bao et al. [7] and Zheng and James [154], who also applied a quasi-static linear elastic model to stress analysis. We use an LCP-based rigid bodysimulator [38] for the purpose of stability and simplicity. We are not discussingrigid body simulation in this paper but refer interested readers to [13] for a com-prehensive review. Figure 4.2 provides a brief illustration of the work flow in ourfracturing system, while detailed major steps in each timestep iteration are summa-rized in Algorithm 3. Each iteration starts with a rigid body simulation, and contactforces are computed for each object. If there are rigid bodies in contact, the contactforces and the gravitational forces will be transformed into the material space, re-spectively, and used as Neumann boundary conditions and inhomogeneous termsfor our boundary integral equation. The dense system generated by such integralequation will be solved using the BiCG-STAB iterative solver [136], combined88Figure 4.2: A graphical visualization of Algorithm 3: our rigid body sim-ulation loop is composed of (a) Rigid Body Dynamics, (b) ResolvingContact, (c) Solving Layer Potential, (d) Stress Analysis and (e) Frac-ture Propagation.with a fast summation method. Once we obtain the layer potential, we can evalu-ate elastic deformation, like displacement and stress, in the material domain. Suchinformation will be used to start and propagate fracture by checking the Rankinecondition. Finally, after all fractures finish propagating or cut completely throughthe material, we update the rigid body list if there are new fragments generated byfracturing process.Due to the characteristics of rigid bodies, linear elasticity is suitable for simu-lating material undergoing small deformations. Even though the quasi-static modeldoesn’t capture elastic wave propagation, which affects fracture as noted by Glonduet al. [42], our indirect boundary integral can be extended to elastodynamic modelsimilarly to acoustic scattering models [58]. Such an extension could be easilycombined with our other modules.894.4 Boundary Element Simulation4.4.1 Elastostatic Constitutive ModelWe assume that the material to be simulated always satisfies hyperelastic andisotropic properties. Traditional approaches involve defining the elastic constitu-tive relation from an energy density and using the first Piola Kirchhoff stress tensor,P. However, for the boundary element method used in this paper, it is necessary tobuild the model upon the second Piola Kirchhoff stress tensor, S. By so doing, weobtain a linearized mapping between the displacement u(x) and the stress tensorS(x), with a small deformation assumption:F(x) = ∇u(x)+ I, x ∈Ωε(x) =12(F(x)+FT (x))− I, x ∈ΩS(x) = C : ε(x), x ∈Ω.(4.1)Here F is the deformation gradient and Ω represents the material space. ε is alinearized right Cauchy-Green strain tensor, and C denotes a symmetric fourthorder material elasticity tensor. The specific linear elasticity model we use isS = 2µε+λ tr(ε)I,where λ , µ are Lame´ parameters.Similar to previous approaches [7, 90, 126], we adopt a quasi-static physicalmodel in the stress analysis, by ignoring the influence of the elastic waves on frac-ture propagation. That is, there are no inertia terms, and we obtain the simplifiedequations∇ ·σ(x) = g, x ∈Ωσ(x) ·n(x) = f(x), x ∈ ∂Ω,(4.2)where g and f represent body forces such as gravity and boundary traction or con-tact force, respectively. Notice that σ represents the Cauchy stress tensor, which is90(a) (b) (c)Figure 4.3: (a) A general material domain enclosed by two codimension oneembedded submanifolds (lines); (b) When the submanifold is exteriorto the enclosed domain, it constructs an interior potential field whosejumping condition on it is positive; (c) Otherwise, it constructs an ex-terior potential field with a negative jumping condition. The surfacenormal is consistent and always pointing out of the material domain.related to S viaσ =1det(F)FSFT . (4.3)Even though in general the stress tensors σ and S are different, it can be shownthat they are equivalent under geometrical linearization when the deformation issmall. By substituting σ into Equation 4.1, we actually obtain the Navier-Cauchyequation, which is solved using a direct boundary integral formulation in [56] toaccomplish interactive haptic rendering.4.4.2 Indirect Boundary IntegralWe now explore the elliptic characteristics of Navier-Cauchy equation and trans-form the PDE into an indirect boundary integral equation using potential theory.A similar idea has been applied to potential flow [23, 64] in computer graphics tosolve Laplace’s equation. In our case, we are dealing with a Poisson-style equationconcerning the inhomogeneous term due to g. We represent the displacement so-lution u(x), when x lies in the material domain Ω enclosed by m codimension oneembedded submanifolds ∂Ωp, as shown in Figure 4.3, as a layer potential with an91additional potential term to account for the body force:u(x) =m∑p=0∮∂ΩpΦ(‖x−y‖2)ρ(y)dy+∮∂Ωpℵ(‖x−y‖2)⊗2 g⊗3 nydy, x ∈Ω(4.4)where the first summand is a single layer potential term while the second is a New-ton potential term [83]. Such a representation only handles Neumann boundaryconditions, as we focus on free-flying rigid body fracturing in this project. Acomplete solution could be achieved by incorporating a double layer potential toaccount for Dirichlet boundary conditions which appear frequently in articulatedbody simulation. The operator⊗n is a conventional n-mode tensor vector multipli-cation, defined as follows using Einstein notation:M ⊗n v =MI1I2···In···IN vInM ∈ RI1×I2×···In···×IN , v ∈ RIn×1(4.5)Φ(r) and ℵ(r) are, respectively, the fundamental solution and a higher order fun-damental solution’s gradient of the Navier-Cauchy equation:Φi j(r) =18piµr[2δi j−(λ +µλ +2µ)(δi j− r,ir, j)]ℵi jk(r) =r,k8piµ− 132piµ(λ +µλ +2µ)χi jk(r)(4.6)whereχi jk(r) = δi jr,k +δikr, j +δ jkr,i+ r,ir, jr,k.δ is the Kronecker delta. ρ(y) denotes the charge density of layer potential inpotential theory and it is the unknown variable in our formulation. To solve forρ(y), we use Equation 4.2 combined with Equation 4.1 and Equation 4.3, which92yieldsf(x) = κ(x)ρ(x)+m∑p=0∮∂Ωp∂ℵ(‖x−y‖2)∂nx⊗2 g⊗3 nydy+12∮∂ΩpC :(Γ(x,y)+ΓT (x,y)) ·nxdy, x ∈ ∂Ωq (4.7)whereκ(x) = sgn(∂Ωq)ω(x)4pi, x ∈ ∂ΩqΓ(x,y) = ∇xΦ(‖x−y‖2)ρ(y), x ∈ ∂Ωq, y ∈ ∂Ωp.The quantity κ(x)ρ(x) accounts for the jump behavior of the second integral oper-ator when x approaches the domain’s boundary; sgn(∂Ω) is a sign function whichequals 1 when ∂Ω encloses an interior domain and −1 otherwise; and ω(x) is asolid angle at the node x. Note that both f(x) and g(x) are represented in the mate-rial space. For a more detailed derivation of the layer potential model we presentedin this section, please refer to Appendix C.4.4.3 Boundary Element DiscretizationIn general the solution of an elliptic PDE can be represented as a single layer poten-tial, double layer potential, or even a mixed potential. Single layer potentials giverise to Fredholm equations of the second kind after taking the normal derivative,which leads to a well-conditioned linear system for Neumann boundary problemafter discretization. Such well-conditioning property is guaranteed by the Fred-holm Alternative theory from functional analysis. In Figure 4.4 (a, b, c, d), wecompare our implementation with another mixed potential method [66] for 2DLaplace’s problem. Both methods achieve similar results for compatible bound-ary conditions; see (a, b). On the other hand, our method generates a smoothersolution in (d), as opposed to (c), for incompatible boundary conditions. We alsoshow that our approach can be easily applied to more complicated domains in (e,f, g, h).We discretize Equation 4.7 using piecewise constant shape functions, and adopt93(a) (b) (c) (d)(e) (f) (g) (h)Figure 4.4: Upper: We compare our formulation and another online imple-mentation by solving the same problems with compatible boundary con-ditions in a and b (ours) as well as with incompatible boundary condi-tion in c and d (ours); Lower: Problems with more complicated do-mains e, g, h with zero Neumann boundary conditions on the circle, 1on one pair of square sides and −1 on the other pair. A solution for asimple square domain f is also included, for comparison.a collocation scheme at triangle centroids:(D+A)ρ = f− g˜ (4.8)where D is diagonal with nonzero entries corresponding to κ(x)ρ(x), and f collectsall the discretized boundary tractions. A is a dense matrix corresponding to the sec-ond summand on the right-hand side of Equation 4.7, while g˜ is composed of thefirst summand. Please refer to Appendix C for assembly details of Equation 4.8.After obtaining the layer potential ρ by solving the above equation, we can evaluatethe deformation displacement using Equation 4.4 and elastic stress by substitutingthe displacement into Equation 4.1 at any point in the material. 3D elastic resultswith displacements evaluated on the boundary surface and volume energy visu-alization are shown in Table 4.1. Other discretization approaches, including the94original mesh deformed mesh elastic energyTable 4.1: Elastic deformation solved by the indirect boundary integral for-mulation. first row: applying a point force load on the top face; secondrow: twisting the cube; third row: stretching the object.Galerkin projection method, Nystro¨m method combined with high order quadra-ture rules like quadrature by expansion (QBX) [63], can also be applied to obtainhigher order convergence rates.During the fracturing process, we need to run stress analysis at candidate nodes,either for the purpose of opening a crack or propagating a crack. We evaluate theelastic stress at these positions and compute the eigendecomposition of it directly,due to the quasi-static elasticity assumption. The maximum positive eigenvalue95is used to decide whether a candidate node should undergo material failure, or itsfracture propagation velocity if it is already a crack node.4.5 Fast Multipole MethodAlthough the sum of D and A in Equation 4.8 forms a well-conditioned matrix, itsdense structure hinders the application of powerful iterative linear system solvers.By making use of a fast summation method like FMM, we can achieve linear run-ning time to solve the dense system; see, for example, the work by Sun et al. [127]for a clear introduction to 2D FMM. The same framework can be extended to the3D version by replacing the Laurent expansion with a spherical multipole expan-sion of the kernel function [44]. Instead of applying such a complicated expan-sion to our second order tensor kernel function, we introduce a simplified kernel-independent FMM which only requires kernel evaluation. Our work is closelyrelated to [153]. Instead of using a layer potential as they did, we make use ofan RBF-style method to approximate far-field interactions. More specifically, forboundary integrals with different integrands (kernel functions), we choose corre-sponding kernel function as the RBF basis for our FMM, so there is no need foreither kernel expansions or extra layer potential evaluations.Given a cluster of particles lying on a smooth manifold, our kernel-independentFMM starts with building an octree spatial partition. For each octree node that hasa non-empty interaction list, we assign two virtual bounding spheres as illustratedin Figure 4.5a. The virtual sphere is identified as either a representing sphere or achecking sphere. The RBF basis is placed on the exterior sphere to approximate theinfluence from far field to local, while the interior sphere is used to approximatethe influence from local to far field. To determine the RBF’s coefficients α , weevaluate the values η on the checking sphere (the first step in Figure 4.5a, b, c, dand e) and solve a linear system Kα = η (the second step in Figure 4.5a, c, d and e).As opposed to the Tikhonov regularization solver used by [153], we observe thatour matrix K is well-conditioned, and therefore we solve the system directly bycomputing the LU decomposition of the matrix. By fixing the size of both virtualspheres relative to the size of the octree node, we pre-compute the LU factorizationfor a cubic box with unit length once at the beginning and apply it with a sizing96(a) Multipole Expan-sion(b) Local Expansion (c) Multipole to Multi-pole(d) Multipole to Local (e) Local to LocalFigure 4.5: A 2D illustration of RBF based (a) multipole expansion (b) localexpansion (c) multipole to multipole translation (d) multipole to localexpansion (e) local to local translation. Circles in blue are checkingcircle while circles in green are representing circle. The red arrowspoint from known data to the unknowns we want to solve. Numbers onthe arrows represent the order of computation.multiplier ϖ during run time. ϖ is the scaling factor of octree node relative to theunit box. We present the RBF approximation of a 2D problem in Figure 4.6.Finally, we follow the conventional FMM scheme and complete our fast sum-mation in three steps. First, we traverse the octree structure in a bottom-up or-der, applying multipole expansions to leaf nodes as well as multipole-to-multipoletranslations to each node (Figure 4.5a, c). The expansion results are stored as RBFcoefficient vectors for each traversed node. Next, a multipole-to-local expansionstep is taken between each tree node and nodes in its interaction list as shown byFigure 4.5d. We finish the summation by traversing top-down through the octree,97(a) Exact Potential (b) RBF Potential (c) Relative ErrorFigure 4.6: An illustration of the use of an RBF method to approximate a2D anisotropic exterior potential field: (a) exact far field generated by1e5 particles with random mass placed nonuniformly inside the centerblack box (b) approximated far field by only 10 degrees of freedomplaced at red dots (representing circle) and red stars are checking circle(c) relative error between exact field and approximated fieldaggregating the far-field influences by applying local-to-local translations to eachnode (Figure 4.5e) as well as an additional composition step to each leaf node. Inthe composition step, we use a local expansion to account for far-field influences(Figure 4.5b) and a direct summation for the influence from near-field particles.Please see Appendix C for a thorough discussion of the implementation details. Aruntime comparison between our FMM and brute force summation is provided; seeTable 4.2 and Figure 4.7.By making use of FMM, we can compute matrix-vector products in linear time.To solve Equation 4.8, we combine FMM with the BiCG-STAB iterative solver andModel Vertices FMM(s) BF(s) Error(L∞)star 2675 5.639 0.095 1.1e-5armadillo 23245 6.522 7.071 2.15e-4fertility 220332 56.642 648.971 7.7e-5moebius 704480 273.339 6577.087 4.3e-5dragon 3609600 949.959 206429.02 6.1e-5Table 4.2: Detailed input information, algorithm timing and accuracy for testsshown in Figure 4.7 run (sequentially) on a desktop with an Intel i7-3770K and 16GB of memory.98Figure 4.7: Runtime comparison between an O(n2) brute force method (red),O(n) complexity (green dash) and our FMM (blue) using the 3D Navier-Cauchy equation’s fundamental solution as an RBF basis.achieve fast convergence as well as an accurate solution. More detailed numericalresults are provided in Section 4.7.4.6 Fracture Propagation TrackingTraditional physically-based fracture simulation generates and evolves cracks bycutting through volumetric mesh elements. Cutting algorithms have been devel-oped for different kinds of elements, including tetrahedra, polyhedra and pointclouds [148]. To avoid volumetric mesh cutting and remeshing, meshless meth-ods have been proposed to track crack evolution in the material space directlywith an additional surface reconstruction or constructive solid geometry (CSG)step [100, 125]. We present a new explicit crack-tracking algorithm that is com-99Figure 4.8: Propagating fracture following a predefined path on sphere (left)and bunny (right). Bottom row shows the mesh quality of crack sur-face.putationally cheaper than volumetric methods, utilizing local triangle mesh oper-ations (see Figure 4.8). Our approach uses a surface tracking library [35], but wehave made significant modifications for crack tracking.4.6.1 Mesh-Based Surface TrackingThe explicit surface tracking approach has become increasingly popular recently,as it preserves more geometric details than implicit methods. Interesting applica-tions include viscoelastic objects [147], water [22] and smoke [23]. We extendthis idea to rigid body brittle fracture simulation by explicitly tracking the crackpropagation in the material space. We make the same assumption as Hegemannet al. [50], namely that material failure starts from the boundary surface. Such anassumption ignores fracturing effects starting in the interior of the material, butis still acceptable for rigid body fracturing led by contact interaction or user cuts,which are the most common scenarios in computer graphics. It would be possiblein future work to extend our method to support interior cracks by adding volumetricevaluation samples and inserting tiny ellipsoid meshes at crack starting positions.1004.6.2 Contact Force Model and Fracture CriteriaAs we use a velocity level LCP-based rigid body solver to resolve rigid body con-tact. The physical outputs are contact impulses Ic instead of contact forces fc. Tocompute the corresponding contact forces, we estimate the contact duration td sim-ilarly to Koschier et al. [68]. Then we assume a constant force acting during thewhole contact period, which leads to a simple contact force model fc = Ic/td . Afterwe transform all contact forces for an object to its material space, we use them astraction boundary conditions for the elastostatic problem of that object and solvefor ρ of the layer potential as described in Section 4.4.3.To open and propagate a crack, we evaluated two different cracking criteria, theRankine condition and Griffith energy minimization. A shape optimization-basedminimal energy approach gives rise to a direct way for crack propagation underthe gradient descent flow, but suffers from inevitable volume loss depending onsurface resolution, which is also evident in Hegemann et al. [50]. By making useof the Rankine condition, we can control the volume loss, which is independent ofsurface resolution, by a user-defined crack thickness. Such parameter can be setvery small to meet the visualization requirement.4.6.3 Fracture Initiation and PropagationAfter solving the layer potential problem, we apply stress analysis to each vertexon the boundary surface. We choose vertices whose principal stress value is largerthan a user-specified threshold, and maximal among all its 1-ring neighbours, ascandidates to open cracks. Once we have a candidate list, we process them sequen-tially. To open a crack at a given vertex, we first compute the cutting line on its1-ring triangle by the crack plane. If the cutting line is almost aligned with a 1-ringedge, we will replace it with the edge in order to avoid a degenerate triangle. Next,we offset the crack plane along its normal and opposite direction by a thicknessparameter specified by the user. We then compute how these two planes intersectthe vertex’s 1-ring edge and apply local remeshing, as illustrated in Figure 4.9:the three red coloring nodes will compose a new crack line, from which the centernode will be an active node which can evolve into the material while the other twowill be inactive at this moment to represent the ends of the crack line. Once a crack101Algorithm 4 Crack Initiation(N )1: clear fracture line listL2: sortN into descending order based on principal stress3: for each node v inN4: if v lies in 1-ring of lines inL5: continue6: else7: l←Open Crack(v)8: add l intoL9: end if10: end for11: returnLline is generated, it will be stored in a crack line list. An overview of our fractureinitiation procedure can be found in Algorithm 4.To evolve a crack line in the material space, we propose a crack velocity modelas follows and use it to propagate each active node on the crack line.vcp = α(h− pro je(h)) (4.9)where h is the unit vector that bisects the angle composed of the correspondingcrack node and its two nearby crack nodes, pro je(·) is the orthogonal projectiononto the principal eigenvector e of stress tensor. The bisection strategy contributesto a relative smooth crack front, as there is no physical rule on how a crack node(a synthetic construct of our geometric representation) should propagate within theplane defined by the principal eigenvector of stress. α is a scalar product of prin-cipal eigenvalue, a control parameter regarding the mesh resolution and a scalingfactor tuned by the user which we set to 0.5. We also check the direction of thepropagation velocity vector to enforce movement into the material. The overlineFigure 4.9: crack initiation first applies local remeshing (A1) and then prop-agates active crack node into the material (A2).102Figure 4.10: crack extension extends the crack line when the end point eithersatisfies Rankine condition or a smoothness criterion.symbol denotes unit vector. As presented in Figure 4.10, the end points of a crackline can become active according to the following criteria: (a) Rankine condition;(b) crack line smoothness. If one end e of the crack line becomes active, we willapply a crack extension operation. We compute the cutting line on the 1-ring tri-angles of e and choose the one which spans a larger angle with the crack line toavoid an inverting crack. Then we apply the same crack plane offset cutting andremeshing operation as described in crack opening operation. Finally, we activatee and use the extended point as the new end of the crack line.As presented in Algorithm 5, a substep time integration will be applied whenwe propagate the crack lines. During each substep, we propagate each crack linesequentially and handle crack line splitting and merging operations. When a crackline is touching the boundary surface, we apply vertex snapping operation similarto [35] and mark the corresponding vertices on the crack line as inactive. Thepropagation process of a crack will stop if all vertices lying on it are inactive. Anintuitive 2D illustration is presented in Figure 4.11. A merge operation occurswhen a crack line forms a loop during the process of propagation as in Figure 4.12.We model this loop using crack circle and apply the same operation to it as crackline in terms of propagation, smoothing and remeshing. If two different crack linesintersect with each other, we apply a passive strategy to deactivate both of them.When a crack circle degenerates to a point or a non-manifold line, we regardthis as material separation and apply a separation operation by a local mesh con-Figure 4.11: crack splitting splits one crack line by deactivating crack nodesthat touch the boundary surface.103Algorithm 5 Crack Propagation(L , sub timestep)1: dt← sub timestep2: for each line l inL3: Update Node Position(l, dt)4: for each active node v in l5: if v hits boundary surface6: label v as inactive7: end if8: end for9: Local Remeshing(l)10: Crack Line Smoothing(l)11: Feature Preserved Smoothing(l)12: Crack Line Extension(l)13: if number of active nodes in l equal 014: Update Mesh Topology(l)15: remove l fromL16: end if17: end fornectivity update (see Figure 4.8 left). We also check if any two nearby nodes ona crack line or a crack circle are labeled as inactive. If so, we apply a topologychange operation to account for crack cutting through material. Given two inactivecrack node, we first use an average crack plane to cut through the material bound-ary surface colored in blue as in Figure 4.13. The crack plane will create a pathon the surface mesh by intersection. A local remeshing operation will be appliedfollowing the intersection path. Next, we open the path and delete the two trianglesFigure 4.12: crack merging turns a crack line into a crack circle when a loopis detected during the propagation.104Figure 4.13: crack topology change is applied to the non-manifold cracknodes. (B1) We use an interpolated cutting plane to locally remesh theboundary surface (blue). (B2) We then delete the triangles sharing theedge of two nearby nodes on the material surface (green) and thickenthe inserted crack on the boundary surface. (B3) Finally, we stitch upthe mesh to make it a closed manifold.sharing the same edge defined by the two inactive crack nodes on the crack surfacecolored in green. Finally, we close up the mesh by connecting the hanging verticeswith new triangles.4.6.4 Remeshing OperationMesh quality is always an essential component of physical simulation. Duringthe substep time integration, a remeshing operation will be applied within eachsubstep, as mentioned in Algorithm 5. Various metrics as well as remeshing ap-proaches have been proposed for volumetric and surface meshes. Our remesh-ing operations are similar to what have been reported in [21, 35], including edgesplitting, edge collapsing, edge flipping. In addition, we make several extra mod-ifications for valid operations to keep the crack line consistent before and afterremeshing, which are listed as follow:• During edge splitting, if one of the edge vertex is non-crack node, the new105vertex is also non-crack node. Otherwise, the new vertex is a crack node.• During edge collapsing, if both edge vertices are crack node or non-cracknode, the merged vertex will be crack node or non-crack node respectively.Otherwise, the merged vertex will be a crack node.• An edge flipping operation will take place if and only if at least one vertexof the candidate edge and one of the new edge are non-crack nodes.Besides the above local remeshing operations, during each step we also apply asmall amount of ridge preserving smoothing and crack line smoothing operations.4.7 Results and DiscussionBoth the single layer potential and Newton potential kernel functions in Equa-tion 4.6 and their gradients are singular. In our implementation, we mollify thesingular kernel by clamping the distance variable r to 10−5m and apply Gaussianquadrature on each triangle when assembling the discretized linear system in Equa-tion 4.8. To verify our physical model, we implement solvers for the 2D Laplacianproblem in Matlab (Figure 4.4) and the 3D elastostatic problem in C++ (Table 4.1).The 2D tearing result, shown in the submission video, is composed of a plane stressmodel implemented in Matlab and a 2D ElTopo tracking library mexed from C++.The linear system is generated by a second kind Fredholm equation and is usuallywell-conditioned. We achieve convergence in less than 50 iterations (mostly 20 to30) using a BiCG-STAB solver for all the examples we show in this paper. Byreplacing the matrix vector multiplication operations with FMM in each iteration,we achieve a linear running time with the model’s surface resolution. In our FMMimplementation, each octree node will be accompanied by two virtual spheres asillustrated in Figure 4.5. For a unit cube centered on the origin, we place bothspheres, interior and exterior, co-centered with the cube. The size of the interiorsphere and exterior sphere are set to 1.2 and 3.0 respectively. This stencil is appliedto each octree node and scaled according to its size.Our fracture propagation method is built on LosTopos, an explicit surface track-ing library. As mentioned in Section 4.6.3 and Section 4.6.4, we modify the106Figure 4.14: Three sculptures on a dinner table are generated using ourmesh based brittle fracturing method. Our explicit tracking approachpreserves tiny but interesting surface details during fracture propaga-tion.remeshing and topological changing strategies to adapt fracture propagation prob-lem. Our current approach doesn’t support adaptivity, so we pre-process all thedownloaded mesh files at the beginning of each simulation (Figure 4.1, Figure 4.15,Figure 4.14, Figure 4.16), which is composed of several steps of remeshing andfeature-preserved smoothing. In Figure 4.8, we show examples of crack propa-gation following a user defined path. The result is rendered in wireframe so themesh updates can be clearly seen when the crack cuts through the model. Duringthe physical simulation, we adopt the crack growth idea that brittle material re-quires significantly higher energy to start a new crack than to lengthen an existingcrack by the same distance, as mentioned in [121], and use two different Rank-ine conditions for fracture initiation and fracture propagation. These two conditionnumbers as well as the crack thickness can be tuned by an artist, even varying themin material space, to achieve the desired fracturing pattern and simulation result.We implemented a velocity level LCP based rigid body solver for results shownFigure 4.15: Crack propagation inside the material due to constant externalforces applied in the opposite direction at the top and bottom face.107in Figure 4.1, Figure 4.14 and Figure 4.16. In Figure 4.14, we apply different Lame´parameters and Rankine condition numbers to each object to achieve a differentfracturing behavior. In Figure 4.1 and Figure 4.16, we demonstrate how the frac-tures propagate in the material space due to the rigid body contacts in the worldspace. We use the explicit Euler method for time integration and adopt a relativelysmall simulation timestep, 0.001s, in order to resolve enough contacts for fracturegeneration. Detailed runtime information of the simulation results is listed in Ta-ble 4.3. Our collision detection method is based on signed distance field versuspoint shell to handle general non-convex shapes. These modules can be replacedby other rigid body dynamic engines based on explicit mesh contact, like Bullet.In the example of Figure 4.15, we apply a pair of constant forces in the oppositedirection at the top and bottom faces of the volumetric shape, and demonstrate howthe crack initiates and propagates inside the material.Limitations and Future Work Since our implementation is very aggressive inaiming to eliminate all volumetric degrees of freedom, there are some limitationsin our proposed method. For example, we don’t consider the case where materialfailure starts in the interior. Secondly, our passive strategy of dealing with crackmerging can be further improved with more understanding of the physical behav-ior of this phenomenon. These limitations are items for future work. Additionalthemes for future investigation include developing more accurate quadrature rulesfor singular kernel integration, mesh adaptivity for explicit surface tracking, de-Case LPSC(s) BiCG-STAB SEC(s) CTC(s)ceramic plate 176.349 8 ∼ 29 3.383 0.877mode I crack 404.492 22 17.008 1.472glass goblet 363.396 4 ∼ 33 8.411 1.92hollow sphere 289.146 4 ∼ 26 3.857 1.043dinner table 397.443 7 ∼ 47 14.736 1.185Table 4.3: Runtime cost details for the 3D fracture simulation cases includedin this paper. LPSC: average layer potential solve cost; BiCG-STAB:BiCG-STAB iterations; SEC: average stress evaluation cost per fracturepropagation step; CTC: average crack tracking cost per fracture propa-gation step.108Figure 4.16: Rigid body brittle fracture simulation of plate (up) and hol-low sphere (down) with respective fracture propagation in the materialspace (second & third column).coupling between the physical simulation’s degrees of freedom and tracking meshresolution, GPU parallization of kernel independent FMM, etc.4.8 ConclusionBy combining an indirect boundary integral formulation, explicit surface track-ing and a kernel-independent fast multipole method, we present the first effectivemethod for rigid body brittle fracture using the boundary surface mesh only, anddemonstrate its merits. Our method provides a novel direction for qualitative frac-ture simulation. It is accurate, and at the same time computationally economical,and it successfully resolves crack evolution in various settings. We have demon-strated the method’s potential, and have pointed out a few directions for furtherexploration and improvement.109Chapter 5Conclusion5.1 SummaryIn this thesis we have explored several computer graphics applications and newalgorithms revolving around optimal mapping. The applications included:• keyframe interpolation for a two-dimensional animation system (with map-pings to establish correspondences between keyframe shapes as well as map-ping from a weighted combination of keyframes to optimal in-between framesof an animation, with or without additional inertial terms for secondaryphysics-based motion),• UV-parameterization (mapping from three-dimensional mesh positions inthe given embedding to a new two-dimensional parameter space),• deformation for geometric modelling (mapping from an undeformed reststate to an optimal low-distortion pose given some constrained boundary po-sitions prescribed by the modeller, in both two and three dimensions),• quasistatic elastic deformation (mapping from an undeformed rest state to aphysically plausible deformation state, according to applied boundary con-ditions, with strain as an output to drive a fracture process).Topology change as part of the problem plays an important role in several of theseapplications, from animating objects which change shape significantly to physics-110based fracture mechanics (and though we did not directly address it, making cutsin a surface to make it topologically equivalent to a disk is often a critical part ofUV-parameterization). In fact, our animation system uses cohesive zone elementsoriginally proposed for fracture mechanics problems where the eventual crack ge-ometry is known in advance, as it is for our animation problems; meanwhile forour fracture simulations we did not assume a known crack path and thus developeda new algorithm based on explicit mesh tracking.The approaches we have taken to solving the optimal mapping problem rangefrom a previously established local-global iteration in Chapter 2, where we focuson the new features of the system (multiple-mesh correspondences, cohesive zonetopology changes, etc.), to a new optimization algorithm for a broad class of hardnonlinear and non-convex problems to drastically improve on local-global iterationin Chapter 3, to the opposite extreme in Chapter 4 where the small displacementcontext meant we could simplify the mapping problem to a linear constant coef-ficient PDE with known fundamental solution, opening up the possibility of anasymptotically faster boundary integral solver.5.2 Future WorkThere are some obvious connections between projects which could be further ex-plored. First and foremost, our motivation for Chapter 3 was in improving thelocal-global iterative solves for Chapter 2, thus checking how beneficial usingBCQN is in the keyframe animation system is natural—although there is a questionas to how to model the cohesive zone elements in the Laplacian.We identified in Chapter 3 that for large 3D problems, direct Cholesky fac-torization of matrices is the major bottleneck since it scales quadratically in run-time cost and superlinearly in storage. We mitigated this both by using a scalarLaplacian (applied to each coordinate independently), with one third the rows andcolumns and one ninth the nonzeros of a Hessian-like matrix used in Newton-related methods, and by using a fixed matrix that only needs factoring once insteadof each iteration. However, the bottleneck remains for BCQN. The Fast MultipoleMethod from Chapter 4 promises O(n logn) solutions to shift-invariant linear ellip-tic problems: we should be able to adapt it to a discrete volumetric summation as111opposed to boundary summation, enabling a much more scalable solve for BCQN.Moreover, the same asymptotics applies to the equations of linear elasticity whichmay prove to be a more effective preconditioner for elasticity-like nonlinear prob-lems than just the Laplacian. On the other hand, it is very difficult to improve on thefloating-point efficiency of modern direct solvers (as a percentage of the theoreticalpeak arithmetic performance of the processor) so it is hard to estimate how largea 3D problem would be necessary for an FMM approach to actually gain in wall-clock performance—though its savings in memory would probably be appreciableon quite small problems, and perhaps even two-dimensional problems.Going the other way, depending on problem size and the desired complexity ofthe deformation mechanics, it may be useful to use BCQN in a volumetric elasticitysolve for fracture, using our mesh-based fracture tracking alongside a geometry-embedded-in-a-coarse-mesh approach to the elasticity to continue to avoid the needfor conforming a 3D tetrahedral mesh to dynamic cracks.Beyond making connections between the projects, there are several more ex-tensions we are interested in. Our current interactive animation editing describedin Chapter 2 relies on energies posed on triangle meshes with linear elements. Thismight be accelerated (and thus open up new possibilities for its use or generaliza-tion) by way of dimensional reduction. Our dynamic interpolation can be viewedas a generalized Lagrangian formulation of elastodynamics, with the elastic po-tential replaced by an averaged potential. According to the Lagrange-d’Alembertprinciple, it can be transformed into the nonlinear partial differential equation sys-temMx¨+Dx˙+∂W I(x,w)∂x= fext ,where M and D are the mass and damping matrix respectively; approximating thiswith a small but well-chosen basis of typical deformations can make time integra-tion much faster.In this thesis, we have mainly focused on problems where map distortion can bemeasured using principal stretches, the singular values of the deformation gradient.There are many more different mapping characterizations that can be utilized, likeBeltrami coefficients (which are popular in holomorphic geometry) or higher orderdifferential operators related to the biharmonic problem. Questions related to the112performance of the BCQN solver may be whether it can perform efficiently in thosesettings, whether the Laplacian operator may need to be replaced by somethingmore closely connected to the energy, and whether further adaptations may benecessary.While our focus in Chapter 3 was purely on the nonlinear minimization withfixed topology and a given initial condition, it is natural to ask whether we couldincorporate topology change (such as the choice of cuts in a surface to make ithomeomorphic to a disc), and whether we can improve on the Tutte embedding forlocally injective initial conditions.113Bibliography[1] N. Aigerman, R. Poranne, and Y. Lipman. Seamless surface mappings.ACM Transactions on Graphics (TOG), 34(4):72, 2015. → pages 48, 49[2] M. Alexa, D. Cohen-Or, and D. Levin. As-rigid-as-possible shapeinterpolation. Proc. ACM SIGGRAPH, 2000. → page 7[3] P. Alliez, D. Cohen-Steiner, M. Yvinec, and M. Desbrun. Variationaltetrahedral meshing. Proc. ACM SIGGRAPH, 24(3):617–625, 2005. →pages 17, 18[4] B. Alpert, G. Beylkin, R. Coifman, and V. Rokhlin. Wavelet-like bases forthe fast solutions of second-kind integral equations. SIAM J. Sci. Comput.,14(1):159–184, Jan. 1993. ISSN 1064-8275. doi:10.1137/0914010. URLhttp://dx.doi.org/10.1137/0914010. → page 87[5] W.-T. Ang. Hypersingular integral equations in fracture analysis. Elsevier,2014. → page 87[6] H. Averbuch-Elor, D. Cohen-Or, and J. Kopf. Smooth image sequences fordata-driven morphing. Computer Graphics Forum, 35(2), 2016. → page 7[7] Z. Bao, J.-M. Hong, J. Teran, and R. Fedkiw. Fracturing rigid materials.IEEE Trans. Vis. Comp. Graph., 13(2):370–378, 2007. → pages85, 86, 88, 90[8] J. Barbicˇ, M. da Silva, and J. Popovic´. Deformable object animation usingreduced optimal control. Proc. ACM SIGGRAPH, 28(3):53:1–53:9, 2009.→ page 9[9] J. Barnes and P. Hut. A hierarchical o(nlogn) force-calculation algorithm.Nature, 324:446–449, 1986. → pages 87, 160[10] W. Baxter, P. Barla, and K. Anjyo. N-way morphing for 2d animation.Comput. Animat. Virtual Worlds, 20, 2009. → page 7114[11] W. Baxter, P. Barla, and K.-i. Anjyo. Compatible embedding for 2d shapeanimation. IEEE Transactions on Visualization and Computer Graphics,15(5):867–879, 2009. → pages 7, 9, 13, 16, 17[12] M. Ben-chen, C. Gotsman, and G. Bunin. Conformal flattening bycurvature prescription and metric scaling. Computer Graphics Forum, 27,2008. → page 48[13] J. Bender, K. Erleben, J. Trinkle, and E. Coumans. Interactive simulationof rigid body dynamics in computer graphics. In Eurographics 2012 Stateof the Art Reports, 2012. → page 88[14] M. Bergou, S. Mathur, M. Wardetzky, and E. Grinspun. Tracks: Towarddirectable thin shells. Proc. ACM SIGGRAPH, 26(3), 2007. → page 33[15] D. P. Bertsekas. Nonlinear Programming. Athena Scientific, 2016. →pages 45, 51, 61[16] G. E. Blandford, A. R. Ingraffea, and J. A. Liggett. Two-dimensional stressintensity factor computations using the boundary element method.International Journal for Numerical Methods in Engineering, 17(3):387–404, 1981. → page 86[17] D. Bommes, M. Campen, H.-C. Ebke, P. Alliez, and L. Kobbelt.Integer-grid maps for reliable quad meshing. ACM Trans. Graph., 32(4):98:1–98:12, July 2013. → page 48[18] J. Bonet and A. Burton. A simple orthotropic, transversely isotropichyperelastic constitutive equation for large strain computations. Computermethods in applied mechanics and engineering, 162(1):151–164, 1998. →page 48[19] S. Bouaziz, S. Martin, T. Liu, L. Kavan, and M. Pauly. Projectivedynamics: Fusing constraint projections for fast simulation. Proc. ACMSIGGRAPH, 33(4):154:1–154:11, 2014. → pages 30, 33, 34[20] A. F. Bower. Applied mechanics of solids. CRC press, 2009. → page 49[21] T. Brochu and R. Bridson. Robust topological operations for dynamicexplicit surfaces. SIAM J. Sci. Comput., 31(4):2472–2493, 2009. → pages33, 85, 87, 105[22] T. Brochu, C. Batty, and R. Bridson. Matching fluid simulation elements tosurface geometry and topology. ACM Trans. Graph. (Proc. SIGGRAPH),29(4):1–9, 2010. → page 100115[23] T. Brochu, T. Keeler, and R. Bridson. Linear-time smoke animation withvortex sheet meshes. In Proc. ACM SIGGRAPH/Eurographics Symp.Comp. Anim., pages 87–95, 2012. → pages 86, 87, 91, 100[24] O. Busaryev, T. K. Dey, and H. Wang. Adaptive fracture simulation ofmulti-layered thin plates. ACM Trans. Graph. (Proc. SIGGRAPH), 32(4):52:1–6, 2013. → page 86[25] I. Chao, U. Pinkall, P. Sanan, and P. Schro¨der. A simple geometric modelfor elastic deformations. ACM Transactions on Graphics (TOG), 29(4):38,2010. → pages 48, 71[26] I. Chao, U. Pinkall, P. Sanan, and P. Schro¨der. A simple geometric modelfor elastic deformations. Proc. ACM SIGGRAPH, 29(4):38:1–38:6, 2010.→ pages 8, 10, 30[27] R. Chen and O. Weber. Gpu-accelerated locally injective shapedeformation. ACM Trans. Graph., 36(6):214:1–214:13, Nov. 2017. →page 50[28] R. Chen, O. Weber, D. Keren, and M. Ben-Chen. Planar shapeinterpolation with bounded distortion. Proc. ACM SIGGRAPH, 32(4):108:1–108:12, 2013. → pages 8, 69[29] Y. Chen, T. A. Davis, W. W. Hager, and S. Rajamanickam. Algorithm 887:Cholmod, supernodal sparse cholesky factorization and update/downdate.ACM Transactions on Mathematical Software (TOMS), 35(3), 2008. →page 68[30] Z. Chen, M. Yao, R. Feng, and H. Wang. Physics-inspired adaptive fracturerefinement. 33(4):113:1–113:7, 2014. → page 86[31] S. Chunrungsikul. Numerical quadrature methods for singular and nearlysingular integrals. PhD thesis, Brunel University, 2001. → page 87[32] S. Claici, M. Bessmeltsev, S. Schaefer, and J. Solomon. Isometry-awarepreconditioning for mesh parameterization. In Proceedings of theSymposium on Geometry Processing, 2017. → page 51[33] E. Corona, P.-G. Martinsson, and D. Zorin. An direct solver for integralequations on the plane. Applied and Computational Harmonic Analysis,38:284 – 317, 2015. → page 169116[34] R. W. Cottle, J.-S. Pang, and R. E. Stone. The Linear ComplementarityProblem. Society for Industrial & Applied Mathematics (SIAM), 2009. →page 61[35] F. Da, C. Batty, and E. Grinspun. Multimaterial mesh-based surfacetracking. ACM Trans. on Graphics (Proc. SIGGRAPH), 33(4):112:1–11,2014. → pages 85, 87, 100, 103, 105[36] B. Dalstein, R. Ronfard, and M. van de Panne. Vector graphics animationwith time-varying topology. Proc. ACM SIGGRAPH, 34(4), 2015. → page7[37] M. Desbrun, M. Meyer, and P. Alliez. Intrinsic parameterizations of surfacemeshes. In Computer Graphics Forum, volume 21, pages 209–218. WileyOnline Library, 2002. → page 48[38] K. Erleben. Velocity-based shock propagation for multibody dynamicsanimation. ACM Trans. Graph., 26(2):12, 2007. → page 88[39] A. Fischer. A special newton-type optimization method. 24:269–284,1992. → page 62[40] A. Frangi. Fracture propagation in 3d by the symmetric galerkin boundaryelement method. International Journal of Fracture, 116(4):313–330, 2002.→ page 86[41] X.-M. Fu, Y. Liu, and B. Guo. Computing locally injective mappings byadvanced mips. Proc. ACM SIGGRAPH, 34(4), 2015. → page 50[42] L. Glondu, M. Marchal, and G. Dumont. Real-time simulation of brittlefracture using modal analysis. IEEE Trans. Vis. Comput. Graph., 19(2):201–209, 2013. → page 89[43] A. Golas, R. Narain, J. Sewall, P. Krajcevski, P. Dubey, and M. Lin.Large-scale fluid simulation using velocity-vorticity domaindecomposition. ACM Trans. Graph., 31(6):148:1–9, 2012. → page 86[44] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J.Comput. Phys., 73(2):325–348, Dec. 1987. ISSN 0021-9991.doi:10.1016/0021-9991(87)90140-9. URLhttp://dx.doi.org/10.1016/0021-9991(87)90140-9. → pages 87, 96, 159[45] M. Gromov. Metric structures for Riemannian and non-Riemannianspaces. Springer Science & Business Media, 2007. → page 7117[46] W. Hackbusch and Z. Nowak. On the fast matrix multiplication in theboundary element method by panel clustering. Numerische Mathematik, 54(4):463–491, 1989. ISSN 0029-599X. doi:10.1007/BF01396324. URLhttp://dx.doi.org/10.1007/BF01396324. → page 87[47] P. Hanrahan, D. Salzman, and L. Aupperle. A rapid hierarchical radiosityalgorithm. In Proc. ACM SIGGRAPH, pages 197–206, 1991. → page 87[48] B. Heeren, M. Rumpf, M. Wardetzky, and B. Wirth. Time-discretegeodesics in the space of shells. Comput. Graph. Forum, 31(5):1755–1764,2012. → pages 7, 8, 10[49] B. Heeren, M. Rumpf, P. Schro¨der, M. Wardetzky, and B. Wirth. Exploringthe geometry of the space of shells. Comput. Graph. Forum, 33(5):247–256, 2014. → page 8[50] J. Hegemann, C. Jiang, C. Schroeder, and J. M. Teran. A level set methodfor ductile fracture. In Proc. ACM SIGGRAPH/Eurographics Symp. Comp.Anim., pages 193–201, 2013. → pages 87, 100, 101, 166[51] N. Higham. Estimating the matrix p-norm. Numerische Mathematik, 62,1992. → page 58[52] R. W. Hockney and J. W. Eastwood. Computer Simulation Using Particles.Taylor & Francis, Inc., Bristol, PA, USA, 1988. ISBN 0-85274-392-0. →page 87[53] K. Hormann and G. Greiner. Mips: An efficient global parametrizationmethod. Technical report, DTIC Document, 2000. → pages 48, 49[54] K. Hormann, G. Greiner, and E.-N. U. G. C. G. GROUP. MIPS: AnEfficient Global Parametrization Method. Defense Technical InformationCenter, 2000. → page 17[55] G. Irving, J. Teran, and R. Fedkiw. Invertible finite elements for robustsimulation of large deformation. In Proc. ACM SIGGRAPH/EurographicsSymp. Comp. Anim., pages 131–140, 2004. → page 85[56] D. L. James and D. K. Pai. ArtDefo: Accurate real time deformableobjects. In Proc. ACM SIGGRAPH, pages 65–72, 1999. → pages 86, 91[57] D. L. James and D. K. Pai. Multiresolution green’s function methods forinteractive simulation of large-scale elastostatic objects. ACM Trans.Graph., 22(1), 2003. → page 87118[58] D. L. James, J. Barbicˇ, and D. K. Pai. Precomputed acoustic transfer:Output-sensitive, accurate sound generation for geometrically complexvibration sources. ACM Trans. Graph. (Proc. SIGGRAPH), 25(3):987–995,2006. → pages 86, 89[59] L. Jiang, R. H. Byrd, E. Eskow, and R. B. Schnabel. A preconditionedl-bfgs algorithm with application to molecular energy minimization.Technical report, DTIC Document, 2004. → page 51[60] C. Kane, E. Repetto, M. Ortiz, and J. Marsden. Finite element analysis ofnonsmooth contact. Comp. Meth. Appl. Mech. Engrg. 180 (1999), 1999. →page 137[61] P. Kaufmann, S. Martin, M. Botsch, E. Grinspun, and M. Gross.Enrichment textures for detailed cutting of shells. ACM Trans. Graph.(Proc. SIGGRAPH), 28(3), 2009. → page 86[62] A. Kaul and J. Rossignac. Solid-interpolating deformations: Constructionand animation of pips. Computer Graphics Forum, 16(1):107 – 115, 1992.→ page 7[63] A. Ko¨ckner, A. Barnett, L. Greengard, and M. O’Neil. Quadrature byexpansion: A new method for the evaluation of layer potentials. J. Comp.Phys., 252:332–349, 2013. → page 95[64] T. Keeler and R. Bridson. Ocean waves animation using boundary integralequations and explicit mesh tracking. In Proc. ACMSIGGRAPH/Eurographics Symp. Comp. Anim., page 9, 2014. → pages86, 87, 91[65] M. Kilian, N. J. Mitra, and H. Pottmann. Geometric modeling in shapespace. Proc. ACM SIGGRAPH, 26(3), 2007. → page 8[66] S. Kirkup and J. Yazdani. A gentle introduction to the Boundary ElementMethod in Matlab/Freemat. In Proc. 10th WSEAS International Conferenceon Mathematical Methods, Computational Techniques and IntelligentSystems, MAMECTIS’08, pages 46–52, 2008. → page 93[67] A. Kort. Computer aided inbetweening. NPAR, 2002. → page 7[68] D. Koschier, S. Lipponer, and J. Bender. Adaptive tetrahedral meshes forbrittle fracture simulation. In Proc. ACM SIGGRAPH/Eurographics Symp.Comp. Anim., 2014. → page 101119[69] S. Z. Kovalsky, N. Aigerman, R. Basri, and Y. Lipman. Large-scalebounded distortion mappings. ACM Trans. Graph, 34(6):191, 2015. →pages 68, 75[70] S. Z. Kovalsky, M. Galun, and Y. Lipman. Accelerated quadratic proxy forgeometric optimization. Proc. ACM SIGGRAPH, 35(4), 2016. → pages50, 51, 53, 56, 68, 69, 73[71] V. Kraevoy and A. Sheffer. Cross-parameterization and compatibleremeshing of 3d models. ACM Transactions on Graphics (TOG), 23(3):861–869, 2004. → page 9[72] Z. Levi and C. Gotsman. Smooth rotation enhanced as-rigid-as-possiblemesh animation. IEEE Transactions on Visualization and ComputerGraphics, 21(2):264–277, 2015. → page 8[73] B. Le´vy, S. Petitjean, N. Ray, and J. Maillot. Least squares conformal mapsfor automatic texture atlas generation. In ACM Transactions on Graphics(TOG), volume 21, pages 362–371. ACM, 2002. → page 48[74] G. Li, L. Yang, S. Wu, W. Tan, X. Chen, and C. Xian. Planar shapeinterpolation using relative velocity fields. Computer Graphics Forum, 37(5), 2013. → page 7[75] Y. Lipman. Bounded distortion mapping spaces for triangular meshes.Proc. ACM SIGGRAPH, 31(4):108:1–108:13, 2012. → pages 9, 13[76] Y. Lipman. Bounded distortion mapping spaces for triangular meshes.ACM Trans. Graph., 31(4):108:1–108:13, July 2012.doi:10.1145/2185520.2185604. → page 48[77] R. Lipton, D. Rose, and R. Targan. Generalized nested dissection. SIAM J.Numer. Anal., 16(2):346–358, 1979. → page 65[78] L. Liu, L. Zhang, Y. Xu, C. Gotsman, and S. J. Gortler. A local/globalapproach to mesh parameterization. In Computer Graphics Forum,volume 27, pages 1495–1504. Wiley Online Library, 2008. → page 48[79] L. Liu, L. Zhang, Y. Xu, C. Gotsman, and S. J. Gortler. A local/globalapproach to mesh parameterization. In Proceedings of the Symposium onGeometry Processing, SGP ’08, 2008. → pages 30, 33[80] T. Liu, S. Bouaziz, and L. Kavan. Quasi-newton methods for real-timesimulation of hyperelastic materials. ACM Transactions on Graphics(TOG), 36(3):23, 2017. → pages 50, 51120[81] S. Martin, B. Thomaszewski, E. Grinspun, and M. Gross. Example-basedelastic materials. Proc. ACM SIGGRAPH, 30(4):72:1–72:8, 2011. → pages8, 11[82] T. Martin, P. Joshi, M. Bergou, and N. Carr. Efficient non-linearoptimization via multi-scale gradient filtering. In Computer GraphicsForum, volume 32, pages 89–100. Wiley Online Library, 2013. → page 50[83] M. Meßner. Time-dependent body forces within a boundary elementformulation. PhD thesis, Technische Universita¨t Graz, 2008. → pages92, 156[84] M. Messner and M. Schanz. An accelerated symmetric time-domainboundary element formulation for elasticity. Engineering Analysis withBoundary Elements, 34(11):944–955, 2010. → page 86[85] N. Molino, Z. Bao, and R. Fedkiw. A virtual node algorithm for changingmesh topology during simulation. ACM Trans. Graph. (Proc. SIGGRAPH),23(3):385–392, 2004. → page 85[86] M. Mooney. A theory of large elastic deformation. Journal of AppliedPhysics, 11(9):582–592, 1940. → page 48[87] P. Mullen, Y. Tong, P. Alliez, and M. Desbrun. Spectral conformalparameterization. In Proceedings of the Symposium on GeometryProcessing, SGP ’08, pages 1487–1494. Eurographics Association, 2008.→ page 48[88] M. Mu¨ller and M. Gross. Interactive virtual materials. In GraphicsInterface 2004, pages 239–246, 2004. → page 85[89] M. Mu¨ller, L. McMillan, J. Dorsey, and R. Jagnow. Real-time simulation ofdeformation and fracture of stiff materials. In Proc. EurographicsWorkshop on Computer Animation and Simulation, pages 113–124, 2001.→ pages 85, 86, 88[90] M. Mu¨ller, N. Chentanez, and T.-Y. Kim. Real time dynamic fracture withvolumetric approximate convex decompositions. ACM Trans. Graph.(Proc. SIGGRAPH), 32(4):115:1–10, 2013. → pages 86, 90[91] Y. Nesterov. A method of solving a convex programming problem withconvergence rate O(1/sqr(k)). Soviet Mathematics Doklady, 27:372–376,1983. → page 50121[92] J. Neuberger. Steepest descent and differential equations. Journal of theMathematical Society of Japan, 37(2):187–195, 1985. → page 50[93] J. Nocedal and S. Wright. Numerical Optimization. Springer Series inOperations Research and Financial Engineering. Springer New York, 2006.→ pages 49, 50, 51, 52, 58, 68, 71[94] A. Norton, G. Turk, B. Bacon, J. Gerth, and P. Sweeney. Animation offracture by physical modeling. The Visual Computer, 7(4):210–219, 1991.→ page 85[95] J. F. O’Brien and J. K. Hodgins. Graphical modeling and animation ofbrittle fracture. In Proc. ACM SIGGRAPH, pages 137–146, 1999. → page85[96] J. F. O’Brien, A. W. Bargteil, and J. K. Hodgins. Graphical modeling andanimation of ductile fracture. In Proc. ACM SIGGRAPH, pages 291–294,2002. → page 85[97] R. Ogden. Large deformation isotropic elasticity – on the correlation oftheory and experiment for incompressible rubberlike solids. In Proceedingsof the Royal Society of London. Series A, Mathematical and PhysicalSciences, pages 565–584, 1972. → page 48[98] S. I. Park and M. J. Kim. Vortex fluid for gaseous phenomena. In Proc.ACM SIGGRAPH/Eurographics Symp. Comp. Anim., pages 261–270,2005. → page 86[99] E. G. Parker and J. F. O’Brien. Real-time deformation and fracture in agame environment. In Proc. ACM SIGGRAPH/Eurographics Symp. Comp.Anim., pages 156–166, 2009. → page 85[100] M. Pauly, R. Keiser, B. Adams, P. Dutre´, M. Gross, and L. J. Guibas.Meshless animation of fracturing solids. ACM Trans. Graph. (Proc.SIGGRAPH), 24(3):957–964, 2005. → pages 85, 99[101] C. G. Petra, O. Schenk, and M. Anitescu. Real-time stochastic optimizationof complex energy systems on high-performance computers. IEEEComputing in Science & Engineering, 16(5), 2014. → page 68[102] C. G. Petra, O. Schenk, M. Lubin, and K. Ga¨rtner. An augmentedincomplete factorization approach for computing the schur complement instochastic optimization. SIAM Journal on Scientific Computing, 36(2),2014. → page 68122[103] T. Pfaff, R. Narain, J. M. de Joya, and J. F. O’Brien. Adaptive tearing andcracking of thin sheets. ACM Trans. Graph. (Proc. SIGGRAPH), 33(4):110:1–9, 2014. → page 86[104] R. Poranne and Y. Lipman. Provably good planar mappings. Proc. ACMSIGGRAPH, 33(4):76:1–76:11, 2014. → page 9[105] A. Portela, M. Aliabadi, and D. Rooke. The dual boundary elementmethod: effective implementation for crack problems. InternationalJournal for Numerical Methods in Engineering, 33(6):1269–1287, 1992.→ page 86[106] E. Praun and H. Hoppe. Spherical parametrization and remeshing. In ACMSIGGRAPH 2003 Papers, SIGGRAPH ’03, pages 340–349, 2003. → page9[107] E. Praun, W. Sweldens, and P. Schro¨der. Consistent meshparameterizations. In Proceedings of the 28th Annual Conference onComputer Graphics and Interactive Techniques, SIGGRAPH ’01, pages179–184, 2001. → page 9[108] M. Rabinovich, R. Poranne, D. Panozzo, and O. Sorkine-Hornung.Scalable locally injective mappings. ACM Trans. Graph., 36(2):16:1–16:16, Apr. 2017. → pages 51, 53[109] M. Rezayat, D. Shippy, and F. Rizzo. On time-harmonic elastic-waveanalysis by the boundary element method for moderate to high frequencies.Computer Methods in Applied Mechanics and Engineering, 55(3):349 –367, 1986. → page 86[110] R. S. Rivlin. Some applications of elasticity theory to rubber engineering.In Proc. Rubber Technology Conference, pages 1–8, 1948. → page 48[111] M. Rumpf and B. Wirth. A nonlinear elastic shape averaging approach.SIAM J. Img. Sci., 2(3):800–833, 2009. → pages 6, 8, 10[112] J. Schreiner, A. Asirvatham, E. Praun, and H. Hoppe. Inter-surfacemapping. In ACM SIGGRAPH 2004 Papers, SIGGRAPH ’04, pages870–877, 2004. → page 9[113] C. Schu¨ller, L. Kavan, D. Panozzo, and O. Sorkine-Hornung. Locallyinjective mappings. Computer Graphics Forum, 32(5):125–135, 2013. →pages 9, 13, 14, 16, 132123[114] J. Shewchuk. Triangle, 2005. URLhttp://www.cs.cmu.edu/∼quake/triangle.html. → pages 13, 131[115] A. Shtengel, R. Poranne, O. Sorkine-Hornung, S. Z. Kovalsky, andY. Lipman. Geometric optimization via composite majorization. ACMTransactions on Graphics (TOG), 36(4):11, 2017. → pages49, 50, 51, 53, 68, 69, 71[116] E. Sifakis, K. G. Der, and R. Fedkiw. Arbitrary cutting of deformabletetrahedralized objects. In Proc. ACM SIGGRAPH/Eurographics Symp.Comp. Anim., pages 73–80, 2007. → page 85[117] J. Sladek, V. Sladek, and S. Atluri. Meshless local petrov-galerkin methodin anisotropic elasticity. Comput Model Eng Sci, 6:477–489, 2004. → page86[118] B. Smith, D. M. Kaufman, E. Vouga, R. Tamstorf, and E. Grinspun.Reflections on simultaneous impact. ACM Transactions on Graphics(Proceedings of SIGGRAPH 2012), 31(4):106:1–106:12, 2012. → page 60[119] J. Smith and S. Schaefer. Bijective parameterization with free boundaries.Proc. ACM SIGGRAPH, 34(4), 2015. → page 18[120] J. Smith and S. Schaefer. Bijective parameterization with free boundaries.ACM Transactions on Graphics (TOG), 34(4):70, 2015. → pages48, 49, 52, 59, 68, 69[121] J. Smith, A. Witkin, and D. Baraff. Fast and controllable simulation of theshattering of brittle objects. Computer Graphics Forum, 20(2):81–90, June2001. → page 107[122] J. Solomon, F. de Goes, G. Peyre´, M. Cuturi, A. Butscher, A. Nguyen,T. Du, and L. Guibas. Convolutional wasserstein distances: Efficientoptimal transportation on geometric domains. Proc. ACM SIGGRAPH, 34(4):66:1–66:11, 2015. → page 7[123] O. Sorkine and M. Alexa. As-rigid-as-possible surface modeling. InSymposium on Geometry processing, volume 4, 2007. → pages 48, 50[124] O. Sorkine and M. Alexa. As-rigid-as-possible surface modeling. InProceedings of the Fifth Eurographics Symposium on GeometryProcessing, SGP ’07, 2007. → pages 18, 30124[125] D. Steinemann, M. A. Otaduy, and M. Gross. Fast arbitrary splitting ofdeforming objects. In Proc. ACM SIGGRAPH/Eurographics Symp. Comp.Anim., pages 63–72, 2006. → pages 85, 99[126] J. Su, C. Schroeder, and R. Fedkiw. Energy stability and fracture for framerate rigid body simulations. In Proc. ACM SIGGRAPH/EurographicsSymp. Comp. Anim., pages 155–164, 2009. → pages 86, 90[127] T. Sun, P. Thamjaroenporn, and C. Zheng. Fast multipole representation ofdiffusion curves and points. ACM Trans. Graph. (Proc. SIGGRAPH 2014),33(4), 2014. → pages 86, 96[128] T. Surazhsky, V. Surazhsky, G. Barequet, and A. Tal. Blending polygonalshapes with different topologies. Computers and Graphics, 25(1):29–39,2001. → page 9[129] V. Surazhsky and C. Gotsman. Controllable morphing of compatible planartriangulations. ACM Transactions on Graphics, 20(4):203–231, 2001.[130] V. Surazhsky and C. Gotsman. Morphing stick figures using optimizedcompatible triangulations. In Proceedings of Pacific Graphics, pages40–49, 2001.[131] V. Surazhsky and C. Gotsman. Intrinsic morphing of compatibletriangulations. International Journal of Shape Modeling, 9(2):191–201,2003.[132] V. Surazhsky and C. Gotsman. High quality compatible triangulations.Eng. with Comput., 20(2):147–156, 2004. → pages 9, 13, 16, 17[133] D. Sy´kora, D. Sedlacek, S. Jinchao, J. Dingliana, and S. Collins. Addingdepth to cartoons using sparse depth (in)equalities. Computer GraphicsForum, 29(2), 2010. → page 7[134] J. Teran, E. Sifakis, G. Irving, and R. Fedkiw. Robust quasistatic finiteelements and flesh simulation. In Proceedings of the 2005 ACMSIGGRAPH/Eurographics symposium on Computer animation, pages181–190. ACM, 2005. → pages 50, 71[135] D. Terzopoulos and K. Fleischer. Modeling inelastic deformation:viscolelasticity, plasticity, fracture. In Proc. ACM SIGGRAPH, volume 22,pages 269–278, 1988. → page 85125[136] H. A. van der Vorst. Bi-cgstab: A fast and smoothly converging variant ofbi-cg for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat.Comput., 13(2):631–644, 1992. → page 88[137] C. Von-Tycowicz, C. Schulz, H.-P. Seidel, and K. Hildebrandt. Real-timenonlinear shape interpolation. ACM Trans. Graph., 34(3):34:1–34:10,2015. → pages 6, 8, 10[138] H. Wang and Y. Yang. Descent methods for elastic body simulation on thegpu. ACM Transactions on Graphics (TOG), 35(6):212, 2016. → page 50[139] O. Weber and D. Zorin. Locally injective parametrization with arbitraryfixed boundaries. Proc. ACM SIGGRAPH, 33(4):75:1–75:12, 2014. →pages 9, 13, 16[140] O. Weber, A. Myles, and D. Zorin. Computing extremal quasiconformalmaps. Computer Graphics Forum, 31(5):1679–1689, 2012. → page 48[141] S. Weissmann and U. Pinkall. Filament-based smoke with vortex sheddingand variational reconnection. ACM Trans. Graph. (Proc. SIGGRAPH), 29(3):115:1–12, 2010. → page 86[142] B. Whited, G. Noris, M. Simmons, R. Sumner, M. Gross, and J. Rossignac.Betweenit: An interactive tool for tight inbetweening. Computer GraphicsForum, 29(2), 2010. → page 7[143] M. Wicke, D. Ritchie, B. M. Klingner, S. Burke, J. R. Shewchuk, and J. F.O’Brien. Dynamic local remeshing for elastoplastic simulation. ACMTrans. Graph. (Proc. SIGGRAPH), 29(4):49:1–11, 2010. → page 85[144] T. Winkler, J. Drieseberg, M. Alexa, and K. Hormann. Multi-scalegeometry interpolation. Computer Graphics Forum, 29(2):309–318, 2010.→ page 8[145] B. Wirth, L. Bar, M. Rumpf, and G. Sapiro. Geodesics in shape space viavariational time discretization. In Energy Minimization Methods inComputer Vision and Pattern Recognition, volume 5681, pages 288–302.2009. → pages 8, 10[146] B. Wirth, L. Bar, M. Rumpf, and G. Sapiro. A continuum mechanicalapproach to geodesics in shape space. International Journal of ComputerVision, 93(3):293–318, 2011. → pages 7, 8126[147] C. Wojtan, N. Thu¨rey, M. Gross, and G. Turk. Deforming meshes that splitand merge. ACM Trans. Graph. (Proc. SIGGRAPH), 28(3):76:1–76:10,2009. → page 100[148] J. Wu, R. Westermann, and C. Dick. Physically-based simulation of cuts indeformable bodies: A survey. In Eurographics 2014 State-of-the-ArtReport, pages 1–19, 2014. → pages 87, 99[149] J. Xing, L.-Y. Wei, T. Shiratori, and K. Yatani. Autocomplete hand-drawnanimations. Proc. ACM SIGGRAPH Asia, 34(6), 2015. → page 7[150] H. Xu, F. Sin, Y. Zhu, and J. Barbicˇ. Nonlinear material design usingprincipal stretches. ACM Transactions on Graphics (TOG), 34(4):75, 2015.→ page 48[151] X. Xu, L. Wan, X. Liu, T.-T. Wong, L. Wang, and C.-S. Leung. Animatinganimal motion from still. Proc. ACM SIGGRAPH Asia, 27(5), 2008. →page 7[152] W. Yang, J. Feng, and X. Wang. Structure preserving manipulation andinterpolation for multi-element 2d shapes. Computer Graphics Forum, 31(7pt2), 2012. → page 7[153] L. Ying, G. Biros, and D. Zorin. A kernel-independent adaptive fastmultipole algorithm in two and three dimensions. J. Comput. Phys., 196(2):591–626, May 2004. ISSN 0021-9991. doi:10.1016/j.jcp.2003.11.021.URL http://dx.doi.org/10.1016/j.jcp.2003.11.021. → pages 87, 96, 162[154] C. Zheng and D. L. James. Rigid-body fracture sound with precomputedsoundbanks. ACM Trans. Graph. (Proc. SIGGRAPH), 29(3), 2010. →pages 85, 86, 87, 88[155] Q. Zhou, J. Panetta, and D. Zorin. Worst-case structural analysis. ACMTrans. Graph., pages 137:1–137:12, 2013. → page 167[156] Y. Zhu, R. Bridson, and C. Greif. Simulating rigid body fracture withsurface meshes. ACM Trans. Graph., 34(4):150:1–150:11, July 2015. →page vi[157] Y. Zhu, J. Popovic´, R. Bridson, and D. M. Kaufman. Planar interpolationwith extreme deformation, topology change and dynamics. ACM Trans.Graph., 36(6):213:1–213:15, Nov. 2017. → page vi127[158] Y. Zhu, R. Bridson, and D. M. Kaufman. Blended cured quasi-newton fordistortion optimization. ACM Trans. Graph., 37(4):40:1–40:14, July 2018.→ page vi128Appendix APlanar Interpolation withExtreme Deformation, TopologyChange and DynamicsA.1 Constructing Multimesh StructureGiven compact non-self-intersecting planar shapes, there are a pair of ambiguitiescharacterized by topological (in-)equivalence. When shapes are topologically in-equivalent, one shape can not be continuously morphed to the other as there is nobi-continuous map (homeomorphism) between them. Thus we cannot build an iso-morphism between the shapes’ standard topology. Nevertheless, one can developa quotient set, P/ ∼, by introducing an equivalence relation, defined as homeo-morphism, to the set of compact non-self-intersecting planar shapes, P. P/ ∼ canthen be generated by two elements, equivalence class of disk, s (simply connected),and annulus, m (multiply connected) under the operations, disjoint union (⊔) andconnected sum (#) as shown in Figure A.1. Thus we require two topological op-erations, cutting and opening, to convert input shapes into an equivalence class inP/∼. Locations and shapes of these cuts and openings are defined on the domain,which then determine both how and where shapes will open and split.When shapes are topologically equivalent, one shape can be smoothly mor-129.../~...[ ][ ][ ][ ][ ] = #Connected Sum[ ] [ ]=ПDisjoint Union[ ] [ ] [ ]Figure A.1: Non-self intersecting planar shapes can be categorized by defin-ing equivalent relationship as homeomorphism. Two shapes are in thesame equivalence class if one can be smoothly morphed to the other.Different equivalence class elements can be converted to the same ele-ment under topological operations, connected sum and disjoint union.phed into another. However, there then exists infinitely many bi-continuous mapsbetween them as both correspondences and path can vary. We provide interactivecontrol of specifying boundary correspondences and inbetweening energy param-eters to allow additional artistic control.Equivalence Tool Our interactive equivalence tool allows users to convert inputartworks to the same equivalence class by outlining and placing desired cuts andopenings. We represent outlines, cuts and openings with a planar straight linegraph. Outlines define the boundary of each drawing’s embedding shape and areequivalent to a planar circle (S1) dividing R2 into two connected components -the bounded inside and the unbounded outside. By specifying orientation of theoutline, clockwise or counter-clockwise, users define which connected componentembeds the drawings.As the specified cuts and openings break the simple closed curve property ofour outlines, we need to reconstruct boundaries after the outline & topo operation(see Figure A.2). We treat the resulting soup of outlines, cuts and openings as a di-rected graph. Orientations of outlines determine directions of their line segments,130Input Drawing Outline & Topo OperationCutOpenbadcefghabcdegfhCorrespondenceBoundary & Cohesive ZoneCohesive Zone (Edge)Compatible MeshingComesh OptimizationFigure A.2: To generate consistent multimesh structure for given input draw-ings, we start with specifying outlines, cuts and openings; then con-struct S1 equivalent boundaries using directed graph flow algorithm;next mark sparse correspondences; and compatibly mesh all shapes todeliver an initial multimesh solution which will be further improved bycomesh optimization.while edges generated by cuts and openings are bi-directed. The final boundariesare reconstructed by applying a directed graph flow algorithm whose details areprovided in Section A.2. After boundary reconstruction, each line segment belong-ing to a cut or opening turns into two duplicated edges with opposing directionsas shown in Figure A.5. We designate each pair of duplicated edges as a cohesiveedge along which we will assign cohesive energies to help form our inbetweeningenergy. We then call collections of cohesive edges a cohesive zone.Boundary Tool Next our boundary tool lets users interactively mark correspon-dences on desired boundaries. Collections of corresponding boundaries acrossshapes are then compatibly reparameterized to obtain pairwise bijectively mappedboundaries suitable for our initial meshing step.Initial Meshing We now have equivalent shapes with desired cuts, openings andcorrespondences specified. Next, we augment our shapes with interior meshes.Given multiple compatible boundaries of planar shapes, we use Triangle [114] to131generate a conforming Delaunay triangulation of one shape. As this meshing stepmay insert Steiner points along the boundaries, we then upsample boundaries of theremaining shapes and ensure that splits in cohesive edges remain consistent. Then,we apply Schueller et al.’s LIM algorithm [113] to map our first shape’s trianglemesh to the remaining example shapes. This will build the initial multimesh struc-ture with mesh topology O = {V,E,T}, cohesive edge sets Ci, boundary vertexsets Bi and mesh geometry Xi (see Figure A.3 for an illustrative example.). As wewill discretize and minimize our inbetweening energy on these meshes, compatiblemeshes with high quality elements and low distortion mappings between shapes arerequired. Thus we apply Comesh Optimization to obtain the final reliable meshes,see Figure A.4.A.2 Directed Graph Flow ImplementationOur interactive equivalence tool allows users to convert example artworks to thesame equivalence class by drawing desired cuts and openings. After an editingsession, we apply a directed graph flow algorithm to construct final boundaries, asthe specified cuts and openings break the simple closed curve property of outlines.We treat the resulting soup of outlines, cuts and openings as a directed graph G .Orientations of outlines determine line segment directions, while edges generatedby cuts and openings are bi-directed.We then apply the directed graph flow algorithm detailed in Algorithm 6. Herewe iteratively look for closed circular flows (directed cycles) until all 1-directed...Figure A.3: An illustrative example of multimesh structure for two shapes:both shape share the same mesh topology but may have different cohe-sive edges, boundary vertice and mesh geometry. Edges colored in redrepresent cohesive edge.132Example B-A-H Example U-SInput:Output:26+4Figure A.4: Comesh Optimization improves both the mapping and meshingquality of the multimesh structure subject to user specified correspon-dences and cohesive edge constraints and keeps consistent mesh topol-ogy throughout the procedure.edges are visited once and all bi-directed edges are visited twice. In detail, hereedge can be visited returns an edge e∗ whose visit counter c(e∗) is not zero,while next edge on graph returns the next edge with respect to e∗ according tothe graph topology and directions field. When there are multiple candidates fornext edge, we choose the one clockwise-closest to e∗. After boundary reconstruc-tion, each line segment belonging to a cut or opening turns into two duplicatededges with opposing directions, see Figure A.5.A.3 Comesh Optimization Energy GradientWe compute the necessary energy gradients here. Gradient of ODT mesh energy is∂D(Xi,O)∂X ji=∑tk∈Γ˜ j(∂ak(Xi)∂X ji·∑vp∈tk ‖X ji −X pi ‖2)8Γ j(Xi), (A.1)133e1e2e3 e4e5e6e7e8e9e10Input Drawinge1e2e3 e4e5e6e7e8e9e10e3 e4e6e7e8e9e10e7e9e10e7+ e7-Outline Cut & OpenDirected Graph Flowb1b2b3Figure A.5: Cuts and openings (red) will break the simple closed curve prop-erty of outlines (green). We apply directed graph flow algorithm toconstruct final boundaries. We treat soup of outlines, cuts and openingsas directed graph, where edges coming from outlines are 1-directed(directed) and edges coming from cuts and openings are bi-directed(undirected). We iteratively look for directed cycles in the graph untilevery 1-directed edges are visited once and bi-directed edges are visitedtwice.and gradient of MIPS map energy is∂C (Xi,X j,O)∂Xqj= ∑tk∈Γ˜q(2∑p=1∂σ∗k∂σ pk∂σ pk (Xi,X j)∂Xqj)ak(Xi),∂C (Xi,X j,O)∂Xqi= ∑tk∈Γ˜q(2∑p=1∂σ∗k∂σ pk∂σ pk (Xi,X j)∂Xqi)ak(Xi)+σ∗k (Xi,X j)∂ak(Xi)∂Xqi,σ∗k =σ1kσ2k+σ2kσ1k(A.2)134where X ji is the j-th vertex in the i-th shape, namely Xi(v j), and Γ˜ j is the one-ring triangles of the j-th vertex. Notice σ∗k (Xi,X j) = σ∗k (X j,Xi) as MIPS energy issymmetric with respect to the source and target shape under the piecewise linearsetting, which can be verified as follows: deformation gradient Fk(Xi,X j) is a linearmap from the tangent space of k-th triangle in shape Xi to tangent space of k-thtriangle in shape X j. For non-degenerate cases, the inverse map Fk(Xi,X j)−1 is thedeformation gradient Fk(X j,Xi) and σpk (Xi,X j) =1σ pk (X j,Xi), p = 1,2.Algorithm 6 Directed Graph Flow1: Inputs:directed graph G = {V,E},E = {e1,e2, · · · ,e f ,e f+1,e f+2, · · · ,eg},i ∈ [1, f ] is 1-directed, i ∈ [ f +1,g] is bi-directed,C = {c(ei)}, i ∈ [1,g].2: Outputs:boundary listB = {bi}.3: Initialize:c(ei) = 1 ∀i ∈ [1, f ],c(ei) = 2 ∀i ∈ [ f +1,g]B = /0,b∗ = /0,e∗ = e1.4: while true do5: if is already in(b∗, e∗)6: add to boundary list(B, b∗)7: b∗ = /08: e∗ = edge can be visited(C )9: if e∗ = /0 break10: else11: add to directed cycle(b∗, e∗)12: c(e∗) = c(e∗)−113: e∗ = next edge on graph(G , e∗)14: end if15: end whileA.4 Continuous Energies for InbetweeningElastic Energy In the continuous setting, we reuse math notation X and x to rep-resent continuous example and deformed shape manifold respectively. Given the135deformation map f : Xi → x from example shape to deformed shape, the ARAPenergy isWD(Xi,x) =∫Xiκ minR∈SO(2)‖d f −R‖2F , (A.3)where df is the corresponding differential map, or the pushforward, R is a rotation,κ is a simple, spatially varying, material stiffness, and ‖ · ‖F is the conventionalFrobenius norm.Cohesive Energy Deformation map f gives a displacement field over Xi whichcan be discontinuous along artist-specified cuts and opens, called cohesive zoneCi ⊂ Xi which is Lebesgue measure 0 subset of Xi, like graphs or nets. Here wereuse the math notation C in the continuous setting. Cohesive energy is designedto measure the gap of discontinuity along Ci,WC(Ci,x) =∫CiE (J x(v) K), v ∈Ci. (A.4)We reuse the math notation v to represent point lying on Ci in the continuous set-ting. The jump operator J·K evaluates the displacement difference between bothsides, x+(v) and x−(v), of cohesive zone,J x(v) K= x+(v)− x−(v), v ∈Ci. (A.5)The energy density E measures the squared norm of this difference to compute theamount of energy accumulated in the zone. Separation occurs when the accumu-lated energy exceeds the separation threshold τ ,E (J x K) =τ˜(J x K), τ˜ < ττ, else. ,τ˜(J x K) = λ ‖J x K‖2, 0≤ λ ,τ.(A.6)To further control cohesive energy, biasing towards greater resistance to shear oralternatively to stretch, we decompose shear and stretch measures of the displace-136ment jump by redefiningτ˜(J x K) = λ [(1−α)‖J x K‖2//+α‖J x K‖2⊥] , α ∈ [0,1]. (A.7)Here ‖ · ‖// and ‖ · ‖⊥ are norms measuring input vector’s tangential and normalcomponent with respect to the direction xˆ of deformed edge pair, defined in Equa-tion A.14, while a change of weights, α , decides their relative contribution.Non-overlap To avoid overlapping during inbetweening, we need to ensure thatpoints x(v) on deformed shape boundaries ∂x do not penetrate the interior. Were-purpose the indicator function from contact mechanics [60]I (x(v)) =0, if x(v)∩ (x\∂x) = /0,∞, else. , v ∈ ∂Xi. (A.8)Here the indicator function I maps boundary point x(v),v ∈ ∂Xi to an extremepenalty energy∞ if it lies in the interior of the deformed shape x and to 0 otherwise.We then integrate up the indicator function over each example shape’s bound-ary to create an energy that ensures boundary points on shape Xi avoid overlap indeformed configuration xWO(∂Xi,x) =∫∂XiI(x(v)). (A.9)A.5 Cohesive Energy DiscretizationIn the discrete setting, cohesive zone is modeled by pairs of overlapping triangleedges, cohesive edge. Since we adopt P1 element to discretize displacement field,displacement jump JxK also varies linearly along cohesive edges. For a cohesiveedge element C ji , we can parameterize displacement jump over it byJx(u)K= (1−u) Jx(v0)K+u Jx(v1)K0≤ u≤1, ∂C ji = {{v+0 ,v−0 },{v+1 ,v−1 }}.(A.10)137As v+ and v− are collocated in example shape Xi, we use v instead to denote them.Given piecewise constant cohesive stiffness λ , the edge element’s cohesive energybelow threshold τ is thenWC(Cji ,x) =∫C jiE (J x(u) K), u ∈ [0,1], τ¯ < τ= λ∫C ji‖J x(u) K‖2= λ∫ l ji0‖J x(u) K‖2dl.(A.11)l ji is the edge length of Cji , which is linearly parameterized asl(u) = u · l ji , 0≤ u≤ 1. (A.12)ThenWC(Cji ,x) = λ∫ l(1)l(0)‖J x(u) K‖2dl(u)= λ∫ 10‖J x(u) K‖2d(u · l ji )= λ l ji∫ 10‖(1−u)Jx(v0)K+uJx(v1)K‖2du= λ l ji {‖Jx(v0)K‖2(13u3−u2+u)+‖Jx(v1)K‖2(13u3)+ Jx(v0)KT Jx(v1)K(−23u3+u2)}∣∣∣∣10=13λ l ji {‖Jx(v0)K‖2+‖Jx(v1)K‖2+ Jx(v0)KT Jx(v1)K}.(A.13)In order to have control over connecting/separating mode behavior, we further splitthe displacement jump JxK into tangential and normal component with respect to xˆ.Even though v+ and v− are collocated in example shape X , it’s not necessarily truefor x(v+) and x(v−) in deformed shape x. For the consistency of math notation con-vention, we use x+(v) and x−(v) instead to represent x(v+) and x(v−) respectively.138Thus xˆ of C ji is defined asxˆ =x+(v0)+ x−(v0)− x+(v1)− x+(v1)‖x+(v0)+ x−(v0)− x+(v1)− x−(v1)‖ (A.14)under the piecewise linear setting. Then we can define ‖ ·‖// and ‖ ·‖⊥ with respectto xˆ as follows,JxK= JxK//+ JxK⊥,JxK// ..= xˆJxKT xˆ, JxK⊥ ..= JxK− xˆJxKT xˆ,‖JxK‖// ..= ‖JxK//‖= (xˆJxKT xˆ)T (xˆJxKT xˆ)= xˆT JxKJxKT xˆ = JxKT (xˆxˆT )JxK,‖JxK‖⊥ ..= ‖JxK⊥‖= JxKT (I2×2− xˆxˆT )JxK.(A.15)By adopting a relative weighting parameter α ∈ [0,1] to control shear/stretch con-tribution and substituting Equation A.15 into Equation A.11, we getWC(Cji ,x) = λ∫C ji(1−α) · ‖Jx(u)K‖2//+α · ‖Jx(u)K‖2⊥. (A.16)Following similar derivations above, we finally achieve the discretization of cohe-sive energy with control over connecting/separating mode behavior.A.6 Transitivity of MIPSWe optimize cyclic summation, instead of full graph summation, of MIPS energyin our comesh algorithm to account for inter-mesh quality. This is due to the factthat MIPS is transitive: if the conformal distortion between Xi and X j, as well asX j and Xk are both small, then the conformal distortion between Xi and Xk is alsosmall. Before we provide the proof for it, we first prove a lemma that will be usedlater.Lemma 1. For matrix A,B ∈ Rn×n, the following inequality holds,‖AB‖F ≤ ‖A‖F‖B‖F . (A.17)139Proof. According to the definition of matrix A’s Frobenius norm, we have‖A‖2F =∑i, j|ai, j|2 =∑i‖Ai∗‖22 =∑j‖A∗ j‖22, (A.18)where Ai∗ is the i-th row vector and A∗ j is the j-th column vector. Due to Cauchy-Schwarz inequality, we have‖Ax‖22 =∑i|Ai∗x|2 ≤∑i‖Ai∗‖22‖x‖22 = ‖A‖2F‖x‖22, (A.19)where x is a column vector. Thus we have‖AB‖2F =∑j‖AB∗ j‖22 ≤∑j‖A‖2F‖B∗ j‖22=‖A‖2F∑j‖B∗ j‖22 = ‖A‖2F‖B‖2F .(A.20)By taking square root of the inequality’s both sides, we prove the lemma.Next we show the transitivity of MIPS, which is defined asσ∗u (Xi,X j) =σ1uσ2u+σ2uσ1u=‖Fu(Xi,X j)‖F‖F−1u (Xi,X j)‖F ,(A.21)where Fu(Xi,X j) is the deformation gradient of triangle tu by deforming from Xi toX j.Lemma 2. If the MIPS distortion between Xi and X j as well as X j and Xk are bothsmall (bounded above), then the MIPS distortion between Xi and Xk is also small(bounded above).Proof. We consider the distortion of a single triangle, tu, among three shapes Xi,140X j and Xk as shown in Figure A.6 and assume the following upper bounds hold,σ∗u (Xi,X j) =‖Fu(Xi,X j)‖F‖F−1u (Xi,X j)‖F ≤ p,σ∗u (X j,Xk) =‖Fu(X j,Xk)‖F‖F−1u (X j,Xk)‖F ≤ q.(A.22)As the deformation gradient Fu(Xi,X j) is a linear map from tu in Xi to tu in X j, itcan be defined asFu(Xi,X j) = Au(X j)A−1u (Xi),Au(Xi),Au(X j) ∈ R2×2,(A.23)for non-degenerate triangles. Au(Xi) and Au(X j) are defined asAu(Xi) =[X1i −X0i ,X2i −X0i],Au(X j) =[X1j −X0j ,X2j −X0j].(A.24)Thus we can rewrite Equation A.22 as‖Au(X j)A−1u (Xi)‖F‖Au(Xi)A−1u (X j)‖F ≤ p,‖Au(Xk)A−1u (X j)‖F‖Au(X j)A−1u (Xk)‖F ≤ q.(A.25)Then we have‖Au(X j)A−1u (Xi)‖F‖Au(Xi)A−1u (X j)‖F·‖Au(Xk)A−1u (X j)‖F‖Au(X j)A−1u (Xk)‖F ≤ pq,⇒ (‖Au(Xk)A−1u (X j)‖F‖Au(X j)A−1u (Xi)‖F)·(‖Au(Xi)A−1u (X j)‖F‖Au(X j)A−1u (Xk)‖F)≤ pq.(A.26)According to 1, we have(‖Au(Xk)A−1u (X j)Au(X j)A−1u (Xi)‖F)·(‖Au(Xi)A−1u (X j)Au(X j)A−1u (Xk)‖F)≤ pq,⇒‖Au(Xk)A−1u (Xi)‖F‖Au(Xi)A−1u (Xk)‖F ≤ pq,⇒ σ∗u (Xi,Xk)≤ pq.(A.27)141Figure A.6: MIPS energy has transitivity property. It is enough to optimizethe cyclic summation instead of the full graph summation. We provethis property by considering chain of maps between single triangle tufrom three shapes, Xi, X j and Xk.Thus the MIPS distortion between Xi and Xk is also bounded.A.7 Automated Blending of Dynamic EffectsWe propose a simple blending control method to balance between hitting keyframesand satisfying dynamics. Instead of utilizing external force to control dynamicslike other space-time optimization approach, our blending method makes use ofnumerical stiffness depending on timestep size hs, material stiffness κ and mate-rial density ρ . We automate the blending of dynamics by reparameterizing theblending weight ξ with a user set compliance control ε ,ξ (s) =ε ·h2s ·κρ · (δ (s)+1) , (A.28)which cancels out the influence of timestep size, material stiffness and materialdensity to turn a multi-parameter control into a single-parameter control. δ (s)measures deformation between the previous frame x¯(w(s−1)) and the target inbe-tweened shape without dynamics Φ(A,w(s)),δ (s) =WD(Φ(A,w(s)), x¯(w(s−1))). (A.29)142Figure A.7: It is easy to blend dynamics or secondary motion into the inbe-tweening using our method to create more expressive and vivid anima-tion. In this flying bird example, we also use inhomogeneous materialby making bird body and wing skeleton much stiffer than head, feet andend of wings region.As the difference between the dynamic and non-dynamic inbetween frames grow,δ (s) increases and the inbetweening x¯(w(s)) is driven towards the interpolatedshape without dynamics Φ(A,w(s)).A.8 Animation ResultsFigure A.8: Our inbetweening method uses local-global solver and interpo-lates shapes with extreme deformation, also shown by Chen et al., with-out any artifacts.Our method can generate high quality results for difficult animation problemsconsisting of extreme deformation (see Figure A.8) and dynamics (see Figure A.7).Our system is easy to use and capable of generating large variety of animation. Toverify this claim, we reproduce the benchmark examples from 12 basic principlesof animation, as shown in Figure A.9, using very few input drawings (at most 4).143SlidingWavingJumpingRollingSquashing&StretchingAligatorFigure A.9: We evaluate our animating method by reproducing animationbenchmark tasks from the 12 basic principles of animation. We achievealmost the same animation results using just few input drawings (2 to4).144A.9 Experiment StatisticsIn this section, we provide detailed statistics tables corresponding to the comeshoptimization evaluation discussed in the paper.A.9.1 Relative Weighting of Energies in Comesh OptimizationObjectiveTable A.1: We take the square-to-hexagon example and evaluate the impactof relative contribution of the ODT mesh and the cyclic MIPS mappingenergies in the objective on our Comesh Optimization’s solution by vary-ing relative weighting as (1 - weight) for meshing and weight for mapenergy. The starting map and mesh energy are 36.394 and 0.516 respec-tively. Note that except for this experiment we use an equal weighting(relative weight of 0.5) in all of our examples and evaluations.Weight Ver(#/#) Total(s) Map Mesh0.13 467/358 2.329 2.114 0.4940.21 467/331 1.881 2.095 0.4940.33 467/340 2.034 2.065 0.4940.44 467/332 3.473 2.054 0.4940.5 467/345 2 2.056 0.4940.62 467/337 2.245 2.052 0.4940.7 467/335 1.186 2.055 0.4940.85 467/327 1.263 2.050 0.4950.93 467/360 1.072 2.046 0.495145A.9.2 Mesh ResolutionTable A.2: We evaluate the scaling behavior of comesh optimization solverby refining mesh resolution of square-hexagon example. As the resolu-tion is doubled, we observe almost linear scaling behavior in terms oftiming.Ver(#/#) Inner(#) Outer(#) Total(s)98/98 59 5 0.065278/277 50 7 0.181545/544 50 8 0.411057/1046 131 13 2.1812087/2074 127 13 4.6914331/4296 117 10 12.2438649/8611 209 17 58.01117131/16963 361 18 220.36235386/35114 311 21 838.02573072/71870 382 26 3358.401165107/164638 447 33 13650.724311522/311183 495 34 51253.368146A.9.3 Number of Example ShapesTable A.3: We evaluate the scaling behavior of comesh optimization solverby increasing number of meshes in the walking character example. Asthe number of meshes increases, we also observe almost linear scalingbehavior in terms of timing.Shape(#) Ver(#/#) Inner(#) Outer(#) Total(s)2 1639/1514 81 6 2.5663 1639/1509 419 13 12.5484 1639/1482 482 12 16.7765 1639/1483 554 16 21.1466 1639/1482 197 7 9.8947 1639/1446 471 16 23.5268 1639/1427 529 12 27.9539 1639/1425 825 11 42.55810 1639/1380 1273 24 77.62111 1639/1380 818 14 85.7112 1639/1357 643 13 56.01413 1639/1342 742 12 94.93314 1639/1314 843 18 82.28315 1639/1337 723 11 104.17116 1639/1357 883 16 78.754147A.9.4 Problem DifficultyTable A.4: We evaluate the performance behavior of comesh optimizationsolver by increasing the problem complexity of pentagon-star example.The solver produces consistently good solutions as we increase varia-tions between pentagon and star. We also observe that similar inputshapes require more time to converge.Ver(#/#) Inner(#) Outer(#) Total(s)hard 1 1683/1516 3149 96 82.906hard 2 1683/1546 1391 53 40.188hard 3 1683/1556 479 13 14.508hard 4 1683/1574 484 22 14.347hard 5 1683/1539 205 9 6.444hard 6 1683/1376 233 9 9.46148Appendix BBlended Cured Quasi-Newton forDistortion OptimizationB.1 EquivalenceTheorem 1. For our energy densities W (σ) = f (σ)/g(σ) with f (σ) > 0 andg(σ)→ 0 as σ → 0, x∗ is a stationary point of {E(x) : a(x) ≥ 0} iff it is a lo-cally injective stationary point of the unconstrained energy E(x).Proof. The 1/g(σ) term drives element energies W(Ft(x))→∞ as at(x)→ 0. Sta-tionary points x∗u of unconstrained E are given by ∇E(x∗u) = 0 and must satisfy|a(x∗u)| > 0. The addition of local injectivity then requires a(x∗u) > 0. Stationarypoints x∗c of {E(x) : a(x)≥ 0} are given by the Karush-Kuhn-Tucker (KKT) con-ditions∇E(x∗c)−∇a(x∗c)λ = 0 and 0≤ λ ⊥ a(x∗c)≥ 0. (B.1)(Here λ = (λ1, ..,λm)T ∈ Rm is a Lagrange multiplier and x ⊥ y is the comple-mentarity condition xtyt = 0, ∀t.) All x∗u satisfy (B.1) with λ > 0. For x∗c sat-isfying (B.1) any λt = 0 =⇒ at(x∗c) = 0 =⇒ W(Ft(x∗c))= ∞. Thus we musthave λ > 0 =⇒ a(x∗c) > 0 so that x∗c are locally injective stationary points of theunconstrained energy E(x).149Vertices Triangles Iteration1.9K 3.1K 19 19 233.5K 6.3K 19 17 376.6K 12.5K 19 27 2412.9K 25.0K 19 26 2425.4K 50.0K 20 19 3350.4K 100.0K 20 25 23100.4K 200.0K 20 19 38197.9K 394.6K 20 28 25435.5K 869.2K 20 29 34880.3K 1,758.1K 21 36 321,650.4K 3,297.5K 21 27 403,221.7K 6,438.7K 21 44 366,386.2K 12,765.6K * * *11,969.0K 23,928.4K * * *CM PN AKVFFigure B.1: AKAP comparison to Newton-type methods. Here we com-pare the convergence performance of the AKAP preconditioner to CMand PN in our UV Parameterization Scaling example (Figure 12 in ourmain paper). We use * to indicate out-of-memory failure for matrixfactorization.150Vertices TrianglesBCQN CM AKVF SLIMIteration Timing(s) Iteration Timing(s) Iteration IterationBull 17.9K 34.5K 169 11.82 25 15.51 25 1.00 148 5.92Camel 40.2K 78.1K 412 92.60 61 103.32 67 1.10 177 2.90Dino 24.6K 47.9K 162 30.43 40 31.40 41 1.03 357 8.93Isis 188.1K 374.3K 347 404.63 33 219.33 37 1.12 † †Cow 3.1K 5.8K 78 1.14 27 1.80 25 0.93 143 5.30Horse 20.6K 39.6K 104 12.93 42 24.38 44 1.05 41 0.98Ratio over CMRatio over CMFigure B.2: AKAP and SLIM comparison. Here we compare the conver-gence performance of SLIM and the AKAP preconditioner in our UVParameterization test example with the ISO energy (Figure 15 in ourmain paper). SLIM and AKAP stencils, and so their fill-in, match CM’s(see Figure 8 in the main paper) and thus require the same per-iterationcompute cost and storage for linear solutions as CM. Here we report thenumber of raw iterations as well as the ratio of iterates per method overCM. We use † to indicate when SLIM does not converge.151Appendix CSimulating Rigid Body Fracturewith Surface MeshesC.1 Boundary Element MethodC.1.1 Indirect Boundary Integral FormulationThe indirect boundary integral formulation is inspired by the potential theory whichis usually applied to Laplace’s equation (scalar or vector). Considering the similar-ity between Laplace and the similarly elliptic linear elastostatic problem, we applythis mathematical tool to stress analysis. This forms a crucial ingredient in ourmesh-based fracture evolution algorithm. Here we briefly introduce the potentialtheory used for Laplace’s problem and then describe how we derive our elastostaticphysical model.C.1.2 Potential TheoryTake a closed codimensional one embedded submanifold ∂Ω, like the contour ofa circle in R2 or a sphere in R3. Let Ω ⊆ Rd be the subdomain enclosed by ∂Ω,either interior or exterior, where d is the dimension of our space. A boundary152integral over such a submanifold has the following form:u(x) =∮∂Ωρ(y)Φ(‖x−y‖2)ds, x ∈Ω,y ∈ ∂Ω. (C.1)Here ds is an infinitesimal patch of ∂Ω, Φ(·) is the kernel function defined inΩ, and ρ(·) is the density function defined on the submanifold. If we choose Φ(·)to be a fundamental solution of Laplace’s equation, then u will be harmonic in Ω.To determine ρ(·), we need to take boundary conditions into account, which in ourcase are of Neumann type. Given a flux function f (·) on ∂Ω, we require∂u∂n(y) = f (y), y ∈ ∂Ω. (C.2)However, to solve for ρ(·), we cannot simply replace the left hand side ofEquation C.2 with a normal derivative of the boundary integral in Equation C.1.The reason is even though the boundary integral in Equation C.1 is continuousthroughout Rd , its normal derivative is continuous only in Ω or Ω¯ - ∂Ω and has ajump condition on the submanifold. Such a jump condition can be stated as∂u∂n(x) = α(x)ρ(x)+∮∂Ωρ(y)∂Φ∂n(‖x−y‖2)dsα(x) = sgn(∂Ω)ω(x)Γ(d2 )2pi d2, x,y ∈ ∂Ω,(C.3)where ω(·) is the solid angle, Γ(·) is the Gamma function. α(x) is ω(x)2pi or ω(x)4piwhen d is 2 or 3 respectively. The function sgn(·) is a sign function, which equals+1 if Ω represents an interior domain and −1 otherwise. Notice the normal nshould be consistently pointing out of the subdomain. By combining Equation C.2and Equation C.3, we can solve for ρ(·). Once ρ(·) is computed, we can evaluateu and ∇u at any point in the space using Equation C.1 or applying the gradientoperator to it.What is described above is called a single layer potential and only handlesLaplace’s problem, while we are more interested in a more general Poisson prob-lem including the constant body force, due to gravity in our case. Now we intro-duce a so-called Newton potential which deals with the inhomogeneous term in153Poisson’s equation. Let’s take a look at Equation C.1 again. If we apply Laplace’soperator to both sides, we will get zero on both sides. Add a domain integral onthe right hand side, likeu(x) =∮∂Ωρ(y)Φ(‖x−y‖2)ds+g∫ΩΦ(‖x− z‖2)dvx,z ∈Ω, y ∈ ∂Ω,(C.4)where g is the constant inhomogeneous term. It is easy to verify that after applyingLaplace’s operator to both sides of Equation C.4, g will remain on the right handside:∇2xu(x) =∮∂Ωρ(y)∇2xΦ(‖x−y‖2)ds+g∫Ω∇2xΦ(‖x− z‖2)dv=∮∂Ωρ(y)δ (‖x−y‖2)ds+g∫Ωδ (‖x− z‖2)dv= 0+g = g.Here δ (·) is the Dirac delta function (or distribution). The extra volume integralis named the Newton potential in Equation C.4. To finish the story, we need to takea further step by transforming the volume integral into a boundary integral, byexploiting a higher order fundamental solution Ψ of Laplace’s equation,∇2Ψ(r) =Φ(r).Then the Newton potential can be transformed into a boundary integral:g∫ΩΦ(‖x− z‖2)dv = g∫Ω∇2Ψ(‖x− z‖2)dv= g∮∂Ω〈∇Ψ(‖x−y‖2),n〉ds.(C.5)Therefore the final indirect boundary integral formulation can be written asu(x) =∮∂Ωρ(y)Φ(‖x−y‖2)ds+g∮∂Ω〈∇Ψ(‖x−y‖2),n〉dsx ∈Ω, y ∈ ∂Ω.154C.1.3 Elastostatic ModelIn our case, we apply the same layer potential idea to the elastostatic problem,representing the displacement field asu(x) =∮∂ΩΦ(‖x−y‖2)ρ(y)ds+∮∂Ωℵ(‖x−y‖2)⊗2 g⊗3 nds,(C.6)whereΦ(·) is a fundamental solution of the Navier-Cauchy equation and⊗n is con-ventional n-mode tensor vector multiplication. The above equation can be viewedas an analogy of Equation C.5 in the elastostatic case, withℵ(·) a third order tensorwhich corresponds to ∇Ψ(·). To compute ρ(·), we make use of the jump conditionagain:∂u∂nx(x) = α(x)ρ(x)+∂∂nx∮∂ΩΦ(‖x−y‖2)ρ(y)ds+∂∂nx∮∂Ωℵ(‖x−y‖2)⊗2 g⊗3 nds, x,y ∈ ∂Ω= α(x)ρ(x)+∮∂Ω∇xΦ(‖x−y‖2)ρ(y) ·nxds+∮∂Ω∂ℵ(‖x−y‖2)∂nx⊗2 g⊗3 nds.The boundary condition we have is an external force f(·) applied on the bound-ary surface, which satisfies the following condition,f(x) = σ(x) ·n(x), x ∈ ∂Ω.Considering the following relations,σ(x) = C : ε(x),ε(x) =12(F(x)+FT (x))− I,F(x) = ∇u(x)+ I,∂u∂nx= ∇u ·nx,(C.7)155we propose the new formula by combining the jump condition and the physicalstationary property,f(x) = α(x)ρ(x)+∮∂Ω∂ℵ(‖x−y‖2)∂nx⊗2 g⊗3 nds+12∮∂ΩC : [∇xΦ(‖x−y‖2)ρ(y)+(∇xΦ(‖x−y‖2)ρ(y))T ] ·nxds.(C.8)C.1.4 DiscretizationIn order to solve for the layer potential, we discretize Equation C.8 using piece-wise constant shape functions and adopt a collocation scheme at triangle centroids,which leads to a dense linear system,(D+A)ρ = f− g˜. (C.9)In the above equation, D is a diagonal matrix whose nonzero entries correspond tothe jump condition, α(x)ρ(x), and f collects all the discretized boundary tractions.A is a dense matrix corresponding to the second boundary integral on the right-hand side of Equation C.8, while g˜ corresponds to the first integral (see paper′sSection 4.3). It is trivial to assemble D and f, and we refer interested readers to[83] for detailed explanation of assembling g˜. In the sequel, we only provide theimplementation details for A.The matrix A is composed of k2 3×3 matrix block Ai j, where k is the number oftriangles in the boundary mesh. Each block Ai j represents the interaction betweentriangles Ti and T j. Given the centroid p and surface normal n of Ti, centroidq and surface area w of T j as well as material properties, including the Poissonratio ν and the shear modulus G, we can compute Ai j as follows. The fundamentalsolution for Navier-Cauchy equation (see paper′s Section 4.2) isΦ(r)i j =116piG(1−ν)r [(3−4ν)δi j + r,ir, j] ,156which can be rewritten asΦ(r) =116piG(1−ν)[(3−4ν)I¯+M¯]M¯ =1r∂ r∂qx· ∂ r∂qx∂ r∂qx· ∂ r∂qy∂ r∂qx· ∂ r∂qz∂ r∂qy· ∂ r∂qx∂ r∂qy· ∂ r∂qy∂ r∂qy· ∂ r∂qz∂ r∂qz· ∂ r∂qx∂ r∂qz· ∂ r∂qy∂ r∂qz· ∂ r∂qzI¯ =1rIr =√(qx−px)2+(qy−py)2+(qz−pz)2∇qr =∂ r∂qx∂ r∂qy∂ r∂qz=qx−pxrqy−pyrqz−pzrr = q−p.Next, we compute ∇pΦ. Since Φ is a second order tensor, its gradient willbe third order. For notational clarity, let us present it in standard form, avoidingEinstein notation. By making use of linearity, we apply the gradient operator toeach column of Φ (viewed as a 3×3 matrix). Then we have∇pΦ(r) =116piG(1−ν)[(3−4ν)∇pI¯+∇pM¯].Since the gradient is a third order tensor, it is useful to describe the process byreferring to individual columns. We use the first column, denoted by col1(Φ), as anillustration example. We have used the symbolic package Mathematica to obtain157the values below:∇pcol1(Φ)(r) =(3−4ν)∇pcol1(I¯)+∇pcol1(M¯)16piG(1−ν)∇pcol1(I¯) =1r3rx ry rz0 0 00 0 0∇pcol1(M¯) =1r32rx− 3r3xr2 −3r2xryr2 −3r2xrzr2ry− 3r2xryr2 rx−3r2yrxr2 −3rxryrzr2ry− 3r2xrzr2 −3rxryrzr2 rx−3r2z rxr2.Similarly, one can compute ∇pcol2(Φ) and ∇pcol3(Φ) analytically. Here we onlyprovide their corresponding component:∇pcol2(I¯) =1r3 0 0 0rx ry rz0 0 0∇pcol3(I¯) =1r3 0 0 00 0 0rx ry rz∇pcol2(M¯) =1r3ry− 3r2xryr2 rx−3r2yrxr2 −3rxryrzr2−3r2yrxr2 2ry−3r3yr2 −3r2yrzr2−3rxryrzr2 rz−3r2yrzr2 ry−3r2z ryr2∇pcol3(M¯) =1r3rz− 3r2xrzr2 −3rxryrzr2 rx−3r2z rxr2−3rxryrzr2 rz−3r2yrzr2 ry−3r2z ryr2−3r2z rxr2 −3r2z ryr2 2rz−3r3zr2Finally, we incorporate the constitutive law shown in Equation C.7 to assemble158Ai j:S1 = 2GE1+2Gνtr(E1)1−2ν IE1 =12[∇pcol1(M¯)+∇pcol1(M¯)T]After computing S2 and S3 in the same way, Ai j is constructed asAi j(p,q,n,w) = w(S1n,S2n,S3n). (C.10)We may now directly construct the matrix A by assembling all 3×3 matrixblocks, and then applying an iterative linear system solver for nonsymmetric sys-tems, such as BiCG-Stab or GMRES. After solving for ρ , we can evaluate the dis-placement and stress information by incorporating Equation C.6 and Equation C.7.C.2 Fast Multipole MethodSolving the dense system described in the last section by applying an iterativesolver may be inefficient. Even though the system is well conditioned (and hencethe iteration count is very small), matrix-vector products take O(n2) floating pointoperations. By replacing the standard matrix-vector multiplication in each iterationwith a fast summation method such as fast multipole method (FMM), an asymp-totic linear running time could be achieved to solve the system. We refer interestedreaders to [44] for a thorough introduction of FMM. Matrix-vector multiplication,Ax = b, can be viewed as an n body problem, where x is a collection of particles’masses that generate a gravity potential field. b stores the potential value evaluatedat corresponding particle’s position. Entry of A, Ai j, represents a gravity potentialinteraction between particle i and particle j. The idea can be generalized by replac-ing the gravity potential kernel function with other non-oscillatory decaying kernelfunctions, like the fundamental solution of Laplace’s equation, the electro-staticequation, the elasto-static equation as well as all the boundary integral kernels usedin our paper. Instead of considering Equation C.10, we stick with gravity potentialpoint of view to give a clear presentation of the algorithm.159Figure C.1: n particlesC.2.1 Algorithm IntroductionGiven n massive particles as shown in Figure C.1, we want to evaluate the grav-ity potential generated by these particles. In particular, we’d like to evaluate thegravity potential at the position of these particles. One can definitely sum up allthe contributions from each particle when considering the gravity potential at eachposition. However, the total complexity is O(n2). A trivial way to reduce the com-plexity is by introducing the concept of center of mass. For example, given twowell separated clusters of n particles, evaluating gravity potential generated by onecluster at the other cluster of particles naively will take O(n2). A more economi-cal approach will approximate the evaluation by first computing the center of massfor each cluster, which costs O(n). Next we take O(1) operations to evaluate thegravity potential generated by one center of mass at the other one. Finally, wespread the computed gravity potentials from the center of mass to all its neighbor-ing particles. The total computational cost for the proposed method is O(n), butthe accuracy depends on how well the clusters are separated. The above methodis used in Barnes-Hut [9] hierarchically and the center of mass approximation iscalled monopole expansion. This simple example provides an essential buildingblock of FMM. A complete description includes also the notion of hierarchicalstructure, multipole expansion, and near-far field decomposition.The FMM algorithm starts with a spatial division tree structure, usually quad160tree or octree as shown in Figure C.2. At each tree level, two cells are regarded aswell-separated if their distance is more than h, where h is the cell size in that level.Otherwise, two cells will be viewed as neighbors. One cell in the quad tree willhave at most 8 neighbors while cell in the octree will have 26 neighbors at most.As in Figure C.3, if we take the green cell as an example, cells colored in red areits neighbors. Furthermore, we introduce the idea of interaction cells colored inblue, which are well-separated from the green cell, but not that far. The formaldefinition can be stated as follows. The interaction cells I (·) are a collection ofcells within the same tree level as green cell g, whose parent cells P(·) (in the treestructure) are neighbors N (·), of the g’s parent cell but they themselves are notneighbors of g:I (g) = {a : P(a) ∈N (P(g))∧a /∈N (g)} .In a quad tree, one cell will have at most 27 interaction cells while 189 in anoctree at most. After the spatial division tree is constructed, FMM will perform thesummation in three steps:• In the first step, we build up a far field approximation for each cell in eachtree level, similar to the evaluation of the center of mass in the above toyexample. However, in order to control the accuracy, FMM uses a multipoleFigure C.2: quad tree Figure C.3: well separated cell161Figure C.4: particle to multipoleexpansion which can be realized in many different ways, such as a sphericalharmonic expansion, layer potential or the RBF-style approach as we pro-pose. The multipole expansion computation starts from the tree leaf cellsup to the tree root. For each leaf cell, we apply a multipole expansion or aparticle-to-multipole step, as shown in Figure C.4. As described in Section5 of the paper, we assign two spheres for each cell where we either place thekernel function or evaluate the potential value. This idea is similar to [153],who classified the spheres as a representing sphere or a checking sphere, de-pending on how they are used. In Figure C.4, the red sphere is representingsphere and the blue sphere is checking sphere. We place k particles uni-formly on each sphere and first evaluate the gravity potential at samples onthe blue sphere generated by the 3 massive particles in the cell:φ p2mi = P1ρ i,where ρ i is the collection of particles’ mass. Then we compute the multipoleexpansion coefficient σ evaluated on the red representing sphere,σ i = P−12 φp2mi .162Figure C.5: multipole to multipoleAfter computing the multipole expansion in the leaf cell level, we pop up themultipole expansion coefficients to their upper level; this process is calledmultipole-to-multipole translation, as shown in Figure C.5. Similar to theparticle-to-multipole step, we use the blue sphere as a checking sphere andthe red sphere as a representing sphere. We first evaluate the gravity poten-tial on the parent cell’s checking sphere generated by all its children cells’representing spheres,φm2mi = M1σ i j +M3σ ik , P( j) =P(k) = i.Then multipole expansion coefficient of the parent cell can be computed asσ i = M−12 φm2mi , M2 = M4We follow the above description and construct the multipole to multipoletranslation level by level.• In the second step, we compute the influence between cells falling in thesame interaction list for each tree level, which is called multipole to localexpansion as shown in Figure C.6. For each cell at any tree level, we aggre-163gate the influence from all its interaction cells,φm2li = T1σ j +T3σ k + · · · .Then the local expansion coefficient τ , which represents far field influence,can be computed asτ i = T−12 φm2li , T2 = T4.Similarly to multipole to multipole translation, we apply this step to eachlevel.• After the second step, we have the influence from interaction cells at eachlevel. In the third step, we pop down such influence from the tree root tothe tree leaf level by level as shown in Figure C.7. Similarly, we start fromevaluating the gravity potential from the parent cell’s representing sphereonto its childrens’ checking sphere,φ l2li j = L1τ iφ l2lik = L3τ i.The popped down local expansion coefficients can be computed for eachFigure C.6: multipole to local164Figure C.7: local to localchildren cell asτ i j = L−12 φl2li jτ ik = L−14 φl2likL2 = L4.The new coefficients will be added to the children’s local expansion coef-ficients computed in the second step and popped down to the next level asabove. Once the local expansion is popped down to the leaf level, we applythe local expansion to particles within each leaf cell as shown in Figure C.8,φ l2pi = Q1τ i.Notice φ l2pi already contains influence from cells well separated from cell ibut the influence from neighbor cells and cell i itself is still not accountedfor, so we finally add the near field contribution by brute froce computationfor each leaf cell as shown in Figure C.9.165Figure C.8: local to particle Figure C.9: near field summa-tionC.3 DiscussionIn this section, we provide a brief discussion on potential extension of both bound-ary element method and fast multipole method.C.3.1 On BEMTo evolve fracture inside the material, various criteria can be applied. We usethe Rankine condition which requires the principal eigenvalue of the stress tensor.Once we obtain ρ by solving Equation C.9, one can evaluate stress information an-alytically at any point inside the material by incorporating Equations C.6 and C.7.The derivation is similar to the previous section, so we are not going to restate ithere. Instead, we give a brief description on how the indirect BEM can be com-bined with another fracture criteria, Griffith energy minimization, which is used in[50]. We show how to analytically compute energy gradient which is the essentialbuilding block in shape optimization. We use the following definition for elasticenergy measure,ψ = σ : ε.Since we use a linearized strain and linear material constitution law, the energywill be a function of deformation gradient F, by which we mean there will be166quadratic terms likeF : F, tr(F)I : Fand a linear term F : I. Thus when we compute ∇ψ , we need to evaluate∇(F : F), ∇(tr(F)I : F), ∇(F : I)Notice the second quadratic term can be rewritten as,tr(F)I : F = tr(F)2.Its gradient turns out to be∇(tr(F)I : F) = 2tr(F)∇tr(F) = 2tr(F)∇(F : I).As the derivation of ∇(F : I) is similar to ∇(F : F), we will just focus on the lattercase,∇(F : F) =∂ (F:F)∂x∂ (F:F)∂y∂ (F:F)∂ z∂ (F : F)∂χ= 2∂F∂χ: F = 2tr[∂F∂χFT ].where χ represents a coordinate parameter. As we have the analytical representa-tion of F based on layer potential, we can also compute ∂F∂χ with little additionaleffort. Thus computing the energy gradient is trivial.Interestingly, it is also possible to compute the gradient of F’s singular values,which might be related to work on stress-aware fabrication, such as [155]. By167applying the gradient to the SVD decomposition of F, we haveF = UΣVT∂F∂χ=∂U∂χΣVT +U∂Σ∂χVT +UΣ∂VT∂χUT∂F∂χV = UT∂U∂χ︸ ︷︷ ︸ωˆUΣ+∂Σ∂χ+Σ∂VT∂χV︸ ︷︷ ︸ωˆVT,where both ωˆU and ωˆVT are skew-symmetric matrices. Since diagonal entries ofmatrix product between skew-symmetric matrix and diagonal matrix will be zero,we have∂σi∂χ= uTi∂F∂χvi.C.3.2 On FMMAs described in the previous section, FMM provides an asymptotically linear com-plexity approach for matrix-vector multiplication, Ax = b, based on n body fastsummation, especially when the summation kernel is non-oscillatory and decay-ing. Below we provide a brief description of FMM from an algebraic point ofview, which relies on matrix low rank approximations.Given n particles as well as spatial division tree structure as shown in Fig-ure C.10, we can reorder the particles based on Morton order such that particlesthat are close in the tree structure will stay close in the algebraic order. We adoptthe definition of neighbor cell, interaction cell and well-separated cell shown inFigure C.3. In particular, we regard interaction between particles lying in neighborcells as full rank and those lying in well-separated cells as low rank. We show thesparse full rank (colored in red) patterns of Figure C.10 for tree levels from 2 to 6in Figure C.12. The low rank setting is due to the usage of non-oscillatory decay-ing kernel function. If we take two clusters of particles and build their interactionmatrix using such kind of kernel function, we will see such low rank property byplotting the matrix’s singular value distribution as shown in Figure C.11. If a ma-trix is rank deficient, the best known rank k approximation is by truncated SVD.Such approaches minimize the error between the original matrix and the approxi-168Figure C.10: particles with depth-6 quad treemation in either L2 or Frobenius norm. However, obtaining such optimal approx-imation not only requires O(n2) matrix construction but also the time consumingSVD decomposition. If we take a look at how FMM works again, the three stepalgorithm actually builds up a hierarchical low rank approximation tree structurefor well-separated particles. Such an algebraic tree structure is similar to hierarchi-cal semi-separable (HSS) or hierarchical block separable (HBS) matrix tree [33].Nevertheless, we view neighbor cell interactions as full rank and want to handlethem directly instead of making approximations as HSS or HBS approaches.To see how FMM can be related to such low rank approximation, let’s taketwo leaf cells of particles as an example. Given the mass ρs of one cell Cs and itsinteraction matrix A with cell Ce, we can evaluate the potential φ e at particles in Ceby matrix-vector multiplication, Aρs = φ e. If these two cells are well-separated, theinteraction can be instead represented as follows based on what we have described169in Section C.2.1:σ s = P−12 φp2ms = P−12 P1ρs(particle to multipole)⇓σP(s) = M−12 φm2mP(s) = M−12 M1σ s(multipole to multipole)⇓...⇓τP◦···◦P(e) =T−12 φm2lP◦···◦P(e) = T−12 T1σP◦···◦P(s)(multipole to local)⇓...⇓τe = L−12 φl2le = L−12 L1τP(e)(local to local)⇓φ l2pe = Q1τe(local to particle).If we collect the above expansions together, we have the following relation,A˜ = Q1L−12︸ ︷︷ ︸L∗1L1·︸︷︷︸L∗2· · · ·T−12︸︷︷︸L∗kT1︸︷︷︸T∗1· · ·M−12 M1︸ ︷︷ ︸M∗2P−12 P1︸ ︷︷ ︸M∗1φ e ≈ φ l2pe = A˜ρs, A≈ A˜,(C.11)where k is the depth of the tree structure and A˜ is a low rank approximation asshown in Figure C.13. If we extend this simple two well-separated clusters inter-action example to general n body problem where full rank interaction and hierar-170Figure C.11: rank deficient interaction between well-separated clusterschical structure need to be taken care of, we obtain a hierarchical matrix algebraicexpansion similar to the HSS or HBS structure as shown in Figure C.14.In the above figure, F∗i stands for full rank block matrix which is composed ofith level neighbor cells interactions. More specifically, F∗1’s matrix blocks representbrute force near field interaction in the bottom level of FMM tree while F∗i+1’s (i ≥1) matrix blocks are (k− i+ 1)th level T∗1 matrix as shown in Equation C.11. Byviewing FMM in this way, we can also observe linear time matrix-vector multipli-cation from an algebraic aspect.171(a) level 2 (b) level 3 (c) level 4(d) level 5 (e) level 6Figure C.12: full rank pattern (colored in red) of Figure C.10 from tree level2 to 6C.4 Numerical ResultsHere we provide some numerical results computed using the above methods. InTable C.1, we apply various surface tractions on different surface models and solvefor displacement and elastic energy on the surface only, whereas in Table C.2, weFigure C.13: low rank approximation of A represented in FMM172Figure C.14: algebraic expansion of FMMsolve for displacement on the surface as well as elastic energy inside the material.173punch cube twist cube compress cubestretch cut punch torus push modelTable C.1: Elastic deformation deformation solved by indirect boundary in-tegral formulation. Color on the surface represents elastic energy.174original meshdeformed meshelastic energyy slicex slicez sliceTable C.2: Elastic deformation solved by indirect boundary integral formula-tion. first row: applying point force load on the top face; second row:twisting the cube; third row: stretching the object.175
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Optimal mapping with topology change for animation and geometry problems
Text