Surface Based Fluid Animation using Integral Equations:Simulation and CompressionbyTodd KeelerB.Sc. Physics, University of Alberta, 2004B.Sc. Computational Science (Math), University of Alberta, 2005M.Sc. Applied and Computational Mathematics, Simon Fraser University, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF OF GRADUATE AND POSTDOCTORALSTUDIES(Computer Science)The University Of British Columbia(Vancouver)September 2017c© Todd Keeler, 2017AbstractThis dissertation looks at exploiting the mathematics of vorticity dynamics and po-tential flow using integral equations to reformulate critical parts of fully dynamicfluid animation methods into surface based problems. These reformulations enablemore efficient calculation and data-structures due to the reduction of the simulationdomain to the two dimensional fluid surface, rather than its volume. We also intro-duce a surface compression and real-time playback method for continuous time-dependent iso-surfaces. This compression method further increases the impact ofour highly efficient surface-based simulation methods.iiLay SummaryThis work contains descriptions of computational algorithms for simulating andcompressing virtual representations of fluid phenomena for their use as visual ef-fects in applications such as feature films, computer games and virtual reality. Thetechniques described are specifically characterized by simulating and compressingdata only described on surfaces such as the boundary between air and water onan ocean or the boundary separating the hot smoky air of a smoke plume with thesurrounding cooler clean air.iiiPrefaceThis thesis contains material that has been published in peer-reviewed forums.Chapter Five is comprised of the paper “Linear-time smoke animation with vor-tex sheet meshes” by Brochu, Keeler and Bridson [14]. The explicit mesh track-ing and vortex sheet implementation were contributed by Brochu while the FastMultipole Method implementation and potential flow boundaries based on integralequations were Keeler’s.Chapter Six contains the paper by Keeler and Bridson “Ocean Waves Anima-tion using Boundary Integral Equations and Explicit Mesh Tracking” [42].Chapter Seven contains the to be submitted paper “ Compact Iso-Surface Rep-resentation and Temporal Compression for Fluid Phenomena ” by Keeler and Brid-son.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Fluid Animation . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 A Mathematical Necessity . . . . . . . . . . . . . . . . . . . . . 42 The Vorticity Formulation of Inviscid Incompressible Flow and itsApplication to Fluid Animation . . . . . . . . . . . . . . . . . . . . . 72.1 Vorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . 9v2.3 Vortex Sheets and Filaments . . . . . . . . . . . . . . . . . . . . 123 Boundary Integral Equations . . . . . . . . . . . . . . . . . . . . . . 133.1 Boundaries and Integrals . . . . . . . . . . . . . . . . . . . . . 133.2 Integral Equations . . . . . . . . . . . . . . . . . . . . . . . . . 143.3 Compact Operators and Fredholm Theory . . . . . . . . . . . . . 163.4 Using Integral Equations for Potential Flow . . . . . . . . . . . . 174 The Fast Multipole Method . . . . . . . . . . . . . . . . . . . . . . . 195 Linear-time Smoke Animation with Vortex Sheet Meshes . . . . . . 235.1 Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 245.3 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 265.4 The Vortex Sheet Mesh . . . . . . . . . . . . . . . . . . . . . . . 295.4.1 Circulation Redistribution . . . . . . . . . . . . . . . . . 315.5 Solid Boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . 325.6 Fast Multipole Method . . . . . . . . . . . . . . . . . . . . . . . 365.7 Rendering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405.8.1 Mesh Simplification . . . . . . . . . . . . . . . . . . . . 405.9 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 416 Ocean Waves Animation using Boundary Integral Equations andExplicit Mesh Tracking . . . . . . . . . . . . . . . . . . . . . . . . . 436.1 Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 436.1.1 Related Work . . . . . . . . . . . . . . . . . . . . . . . . 456.2 Deep Water Waves . . . . . . . . . . . . . . . . . . . . . . . . . 466.2.1 Linearization . . . . . . . . . . . . . . . . . . . . . . . . 486.2.2 Boundary Integral Formulation . . . . . . . . . . . . . . . 486.3 Numerical Method for Boundary Integral Waves . . . . . . . . . . 506.3.1 Step 1: Solving our BIEs by Collocation . . . . . . . . . . 50vi6.3.2 Step 2: Computing ∇Φ . . . . . . . . . . . . . . . . . . 546.3.3 Steps 3 and 4: Time Integration . . . . . . . . . . . . . . 556.3.4 Solid Boundaries and Mesh Adaptivity . . . . . . . . . . 566.4 FMM and GPU acceleration . . . . . . . . . . . . . . . . . . . . 566.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 606.6 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 627 Compact Iso-Surface Representation and Temporal Compression forFluid Phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . 637.1 Synopsis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 637.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 647.3 Lagrangian Surface Representation . . . . . . . . . . . . . . . . . 667.3.1 Pointwise Reconstruction . . . . . . . . . . . . . . . . . 667.3.2 Global Distribution and Reconstruction . . . . . . . . . . 677.3.3 Sharp Features . . . . . . . . . . . . . . . . . . . . . . . 697.3.4 Particle Distribution . . . . . . . . . . . . . . . . . . . . 707.3.5 Spatial Compression . . . . . . . . . . . . . . . . . . . . 717.3.6 Lagrangian Transport . . . . . . . . . . . . . . . . . . . . 717.4 Temporal Compression . . . . . . . . . . . . . . . . . . . . . . . 727.5 Decompression and Re-Meshing . . . . . . . . . . . . . . . . . . 737.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 757.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 798 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83viiList of TablesTable 6.1 Seconds per BiCGSTAB iteration using FMM or GPU acceler-ation for mesh dimensions (Dims) and total number of vertices(Verts) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58Table 7.1 Spatial and temporal compression for a moving smoothed cube. 78Table 7.2 Timings for spatial and temporal Compression and temporalDecompression. . . . . . . . . . . . . . . . . . . . . . . . . . 78Table 7.3 Moving smoothed cube per-frame file sizes in bytes comparedto Corto. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78viiiList of FiguresFigure 1.1 Example of computer generated ocean and water from the fea-ture film “In the Heart of the Sea”, c©Warner Bros. . . . . . . 2Figure 5.1 We represent smoke with an adaptive triangle mesh both forlinear-time simulation, as a vortex sheet, and linear-time inter-active rendering as the boundary of the smoky region. Here adeforming torus influences the buoyant smoke flowing aroundand through it. . . . . . . . . . . . . . . . . . . . . . . . . . 25Figure 5.2 Edge split operation. New edges AE and EB are assigned cir-culation from original edge AB. . . . . . . . . . . . . . . . . 32Figure 5.3 Edge merge operation. Left: The red and green edges andtheir incident triangles are removed from the mesh, leavingtwo quad-shaped holes. Right: The resulting holes are zip-pered with new triangles. Suppose the red and green edges inthe right diagram are the closest matches to the correspondingred and green edges in the left diagram (comparing Euclideandistances of associated edge vectors). These new edges arethen assigned the circulation values from the deleted edges. . . 33Figure 5.4 Simulation time per frame for computing velocity via the Biot-Savart law, for both the direct method and Fast Multipole Method. 37ixFigure 5.5 Simulation time per frame for computing the single layer po-tential on a sphere using BiCGSTAB, computing the matrixmultiply with both the direct method and Fast Multipole Method.Note that the direct method is nearly exactly O(N2) showingthe O(1) convergence of BiCGSTAB . . . . . . . . . . . . . 38Figure 5.6 Our proof-of-concept volume renderer produces, from left toright on the top row, an alpha channel from the smoke mesh,a scattering channel with self-shadowing, and an image-spaceage channel. Below on the left is a sharp composite of the firsttwo, and below on the right includes an age-dependent blur,simulating diffusion in image-space at render-time. . . . . . . 39Figure 5.7 Geometry creation when simulating a smoke plume withouttopological changes (Without TC) and with topological changes(With TC). . . . . . . . . . . . . . . . . . . . . . . . . . . . 41Figure 6.1 Deep Ocean Waves . . . . . . . . . . . . . . . . . . . . . . . 43Figure 6.2 The geometry of our wave formulation. . . . . . . . . . . . . 47Figure 6.3 Ocean Waves initialized with Phillips spectrum [71]. . . . . . 52Figure 6.4 Sequence of images showing waves with longer wavelengthovertaking those with smaller, a characteristic of deep waterdispersion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 6.5 Top down view of the computational mesh showing naturaladaptivity induced by Lagrangian motion. . . . . . . . . . . 58Figure 6.6 Logarithmic plot of table 6.1 showing the number of verticesversus time per iteration asymptotics of the FMM and the GPU. 59Figure 6.7 Waves incident upon a near-vertical wall. . . . . . . . . . . . 60Figure 6.8 Moving balls with wake. . . . . . . . . . . . . . . . . . . . . 60Figure 6.9 Comparison of our waves (top) with a Tessendorf/FFT height-field wave simulation (bottom) using the same initial condi-tions. Coloring represents the potential values on the surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61xFigure 7.1 Resolving sharp features on a cube. On the left we are usingour weighted normal averaging. On the right, no sharp featureresolving is applied. . . . . . . . . . . . . . . . . . . . . . . 70Figure 7.2 The meshed torus, white points show the least squares posi-tions and extruding red lines showing the analytic normal di-rection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 7.3 Convergence of the least squares representation on the surfaceof an analytic torus. The log plot shows the error convergencewith a linear fit for the approximate asymptotic order. . . . . 77Figure 7.4 Smoothed cube test geometry. . . . . . . . . . . . . . . . . . 79Figure 7.5 LSQ points representations of a smoothed cube. . . . . . . . 80Figure 7.6 LSQ points representation of a fluid simulation created fromcommercial software. . . . . . . . . . . . . . . . . . . . . . 80xiAcknowledgmentsThis thesis would not have been possible without the patience and inspired direc-tion of my supervisor, Dr. Robert Bridson.I would also like to acknowledge the candor and kind oversight of my supervi-sory committee, Dr. Dinesh Pai and Dr. Brian Wetton.The MITACS organization played a significant role in providing funding andenriching industrial experience throughout the course of this degree, for which Iam also grateful.Ultimately, the greatest support I received was the patience and loving supportof my family, especially my wife, to whom this thesis is dedicated.xiiDedicationTo my wife Lauri and our four sons, Jacob, Liam, Camden, and Nathaniel.xiiiChapter 1Introduction1.1 Fluid AnimationFluid animation is the art of using computer algorithms to create virtual represen-tations of the motion, geometry, and visualization of fluidic phenomena such aswater or smoke. Fluid animation is used to create visual effects for industries suchas film, computer games, and virtual reality. For feature films, a fluid animationtechnique needs to be directable by an artist as to the movement and placement ofthe effect, but also resemble the physical laws of nature. The animation techniqueneeds to approach a level of believability such that the audience is convinced thatthe simulation was an actual physical reality. These algorithms often push the lim-its of computing power and data storage in an attempt to create a highly-detailedrepresentation of fluid phenomena. They are performed offline and the total com-putation time can take hours or days. Increasing the efficiency and controllabilityof these simulations reduces the amount of artist time spent on creating them andthus provides the impetus for creating more efficient and controllable algorithms.In the context of computer games or virtual reality, however, the speed of theeffect becomes more important. The animation technique or the playback of itsresults must be executed in real-time. This requires a per-frame computation timeof tens of milliseconds or less. Even simple playback of the results of standardanimation techniques in the film industry can exceed these computation limits.Ultimately, games and virtual reality require the fluid effect to be interactive with1the player or environment.Given the requirement that a fluid technique appear realistic, it is no surprisethat successful modern fluid animation techniques are often based on techniquesoriginally developed by the scientific community. These scientific techniques aredesigned to compute numerical approximations of the mathematical models de-scribing fluid motion with rigorous requirements on the approximation error andstability of the technique. Replacing these rigorous requirements with the moresubjective ones required by an animation technique provides a wide design spacein which to research new algorithms that are efficient in creating visually appealingdetail and offer more artistic control. Common innovations in this area include theintegration of cheaper procedural elements, adapting the resolution of the compu-tational domain of the simulation to be dictated by viewing angle, and reduction ofmodel complexity and numerical implementation to achieve greater efficiency. Forfilm, a large scale fluid effect will generally combine both physically accurate ani-mation techniques and more artistic and scientifically inaccurate methods to createthe desired effect.Figure 1.1: Example of computer generated ocean and water from the featurefilm “In the Heart of the Sea”, c©Warner Bros.2Figure 1.1 shows a scene from a feature film that incorporates multiple fluidanimation techniques. All of the water in the scene was generated by a computer.The ocean waves extending to the horizon and around the scene objects were cre-ated using a procedural technique attributed to Tessendorf [71]. The more dynamicfluid simulations around the whales and boat were generated using an industrystandard fluid solver incorporating particle motion with a computational grid com-monly known as a FLIP simulation, originally created by Zhu and Bridson [83].The white water spray in the scene is a simple ballistic simulation artistically con-trolled and sourced from the underlying more physically based simulations.1.2 OverviewMany fluid animation techniques are volumetric, requiring a discrete representa-tion of the fluid volume of interest; those that are surface simulations tend to beeither procedural, such as Tessendorf-style FFT waves, or reduced models suchas shallow or linear wave equations. A critical problem with scaling volumetricmethods is the inherent axial O(n3) complexity of filling three dimensional space.This dissertation presents exploiting the mathematics of vorticity dynamics andpotential flow using integral equations to reformulate critical parts of fully dynamicfluid animation methods into surface based problems. These reformulations allowcalculation and data-structures to be reduced to the two dimensional fluid surface,rather then its volume. This dimensional reduction allows us to circumvent someof the spatial complexity of volumetric simulation while still achieving their levelof dynamic interaction and depth of modelling.As with any kind of reduction of algorithmic complexity, the increased effi-ciency of our methods enables the creation of greater amounts of data - in ourcase, surface data. In order to alleviate the bandwidth associated with using thisdata, we introduce a surface compression and real-time playback method for con-tinuous time-dependent iso-surfaces. This technique enables the playback of thehigher resolution achievable by our surface techniques with lower bandwidth anddisk-storage, further enhancing their desirability.The second through fourth chapters of the dissertation give an intuitive overviewof the mathematical and physical underpinnings of our surface based animation3techniques - vortex sheet dynamics, integral equations, and the Fast MultipoleMethod (FMM). These chapters are meant to be a resource for computer scienceresearchers and graduates, identifying some of the confusion in animation publi-cation and noting relevant sources in mathematical literature for those seeking adeeper understanding of the techniques.The final three chapters are the manuscript portion of the dissertation. Here wepresent the three papers substantiating these novel surface based methods:• “Linear time Vortex Sheet Smoke Simulation” [14]• “Surface Water Waves Using Integral Equations” [42]• “Compact Iso-Surface Representation and Temporal Compression for Real-time playback of Fluid Phenomena” ( to be submitted )1.3 A Mathematical NecessityMost of this thesis depends on a certain level of understanding of Partial Differ-ential Equations. This preliminary section contains a small but critical amount oftheory that is used in every chapter but the last. We therefore give a brief non-rigorous overview for those readers lacking this fundamental.We start with the homogeneous Laplace equation in three dimensions,∇2Φ(~x) = 0. (1.1)We seek a solution for the inhomogeneous Laplace equation in open three di-mensional space, i.e. the infinite domain containing all~x ∈R3 with no boundaries.∇2Φ(~x) = δ (~x). (1.2)The δ is the Dirac delta, a function-like distribution that is zero everywhere butat~0, and has the characteristic that∫R3 δ (~x) = 1.The Dirac delta only really makes sense when convolved with another function,4∫R3f (~x)δ (~x−~y)d~x = f (~y). (1.3)The solution to (1.2)Φ(~x) =14pi|~x| , (1.4)is the called the fundamental solution for Laplace’s equation.This solution is interesting because it actually solves the homogeneous equa-tion everywhere but at~x = 0. In physics, this function notably describes the poten-tial of an electrical point charge.We can now leverage this specific solution to (1.2) to represent generalizedsolutions to Poisson’s equation (i.e. the inhomogeneous extension of Laplace’sequation).∇2Φ(~x) = f (~x). (1.5)Convolving both sides of equation (1.5) with the fundamental solution and us-ing integration by parts requiring the domain where f (x)! = 0 is bounded, we find∫R3f (~x)4pi|~x−~y|d~x =∫R3∇2Φ4pi|~x−~y| d~x (1.6)=∫R3∇2(14pi|~x−~y|)Φ d~x =∫R3δ (~x−~y)Φd~x = Φ(~y). (1.7)Put simply, the solution to (1.5) isΦ(~x) =∫R3f (~y)4pi|~y−~x|d~y. (1.8)This function shows that we can produce the solution to Poisson’s equation ata point, with some requirements on f , simply by a convolution over f with thefundamental solution. Again, in physics this function represents the electrostatic5potential due to a distribution of charge f , and its gradient the electric field. In fact,the theory we have just described is often referred to as “Potential Theory” [43].This integral is used in chapter 2 and 5 to derive the Biot-Savart law for vorticity tovelocity calculations. Boundary integrals create the solution to Laplace’s equationby calculating a distribution of charge across boundaries. We present some ofthis theory in chapter 3, then use it in both chapters 5 and 6. The Fast MultipoleMethod is an O(N) hierarchical summation algorithm for summing the discretepoint distributions resulting from applying quadrature to (1.8). We briefly describeit in chapter 4 and use it again in chapters 5 and 6.6Chapter 2The Vorticity Formulation ofInviscid Incompressible Flow andits Application to FluidAnimationThis chapter is meant to be an overview of the mathematical underpinnings of mod-elling inviscid fluid using vorticity. The aim is to clarify some of the subtleties ofvortex methods. We will derive some of the equations usually given as axiomatic,looking specifically to understand boundary conditions.2.1 VorticityVorticity based animation techniques are primarily used to create smoke effects.In this case there is one fluid, air, of constant density ρ0 with the fluid domainboundaries given by the surface of solid objects immersed in the unbounded fluidregion. More precisely, the aim of vorticity based fluid animation is to produce afast and efficient numerical simulation that approximates visually interesting fluidflow in an unbounded domain Ω = R3/∪iΓi with Γi a finite set of non-overlappingbounded regions ofR3 and ∂Ω=∪i∂Γi. This is done by constructing the numericalsolution of the partial differential equations (PDEs) modelling incompressible flow,7the incompressible Euler equations.For fluid velocity ~u, pressure p, density ρ0, and external forces ~F the incom-pressible Euler equations are∂~u∂ t+(~u ·∇)~u = −∇pρ0+~F (2.1)∇ ·~u = 0. (2.2)Solid boundaries require that the incident fluid velocity have zero normal com-ponent. Thus for solid boundary with solid velocity ~usolid, ~usolid ·~n =~u ·~n where~nis the unit normal on ∂Γ. In addition, for the unbounded domain we require that~u→ 0 as ~x→ ∞. In general, the computationally intensive part of fluid anima-tion when directly approximating the incompressible Euler equations numericallyis satisfying the divergence-free condition (2.2). One of the primary strengths ofvorticity modelling is that through some reformulation we can satisfy this conditioninherently.Vorticity is simply defined as curl of the velocity, ~ω =∇×~u. We will considerflow inR3 which approaches zero at infinity and subsequently add a term to accountfor a solid boundary Γ. Taking the curl of the first of the incompressible Eulerequations, and using some vector identities, (see Saffman’s book, [64]) one caneliminate the pressure dependency and arrive at the vorticity equation∂~ω∂ t+(~u ·∇)~ω = (~ω ·∇)~u+∇×~F . (2.3)In R3, a divergence-free velocity field can always be defined as the curl of a vectorvalued stream function ~Ψ,~u = ∇×~Ψ. (2.4)This leads to deriving the stream function in terms of vorticity,~ω = ∇×∇×~Ψ,~ω = ∇(∇ ·~Ψ)−∇2~Ψ. (2.5)8The choice of ψ is not unique, since we can always add to it the gradient of a scalarfunction, which will not affect the velocity in equation (2.4). We can thereforealways modify Ψ by a gradient of a scalar function, and enforce that ∇ ·~Ψ = 0.Doing this eliminates the first right hand term in equation (2.5). We then can takethe stream function as the solution to a vector Poisson equation with a vorticitysource term.−~ω = ∇2~Ψ, Ψ→ 0 as ~x→ ∞ (2.6)Solving the Poisson problem via convolution with the Green’s function forLaplace’s equation and then taking the curl leads to the following equation forconstructing the velocity from a vorticity field~u =∫R3~ω(~x′)× (~x−~x′)4pi|~x−~x′|3 dV′ (2.7)Equation (2.7) is called the Biot-Savart law. Given an initial distribution of vor-ticity, the Bio-Savart law (2.7) and the vorticity equation (2.3) describe that distri-bution’s evolution and the subsequent time-dependent infinite incompressible flowfield. These equations are the general starting point for most vorticity based anima-tion papers. The subtle part of their construction lies in the fact that the Biot-Savartlaw is derived using the free-space Greens function for Laplace’s equation. Thequestion must be then asked, is the Biot-Savart law valid for domains with solidboundaries, as is assumed by and large in the animation literature? The response isa qualified one, and we will examine it in the next section.2.2 Boundary ConditionsWhen we add solid boundaries to our infinite domain, continuous divergence-freevelocity fields are no longer uniquely defined by the curl of a vector potential.For domains which include solid boundaries and remain simply connected, theHelmholtz-Hodge decomposition gives that any continuous vector field ~u can bewritten as the sum of the curl of a vector potential ~Ψ and scalar potential φ .~u = ∇×~Ψ+∇φ (2.8)9This decomposition isn’t unique, though it can be made so under certain bound-ary conditions. We, however, are more interested in requiring Ψ to be exactly thatvector streamfunction defined in the previous section under the same notation. Byrequiring the above equation to be divergence-free, we arrive at the following re-strictions on φ :∇ ·~u = ∇ ·∇×~Ψ+∇ ·∇φ= ∇2φ = 0 (2.9)The primary requirement for φ then, is that it solves the homogeneous Laplaceequation. From now on, we will consider φ as always satisfying this condition. Thevector field defined by∇φ is called potential flow and is given significant treatmentin introductory fluid textbooks, where in two dimensions it can be described usingcomplex analysis.Since ∇×∇φ = 0, one can add any potential flow to a stream-function velocityfield without modifying the resulting vorticity. This null-space allows us to satisfyno-through boundary conditions by adding the potential flow required to cancelany normal flow on the boundaries resulting from the Biot-Savart law. In our free-space representation without boundaries, we require ~u→ 0 as~x→ ∞. Given solidboundary velocity ~usolid and stream-function velocity values ~ustream , we have thefollowing solution for φ∇2φ = 0 (2.10)∂φ(~x)∂~n=−(~usolid+~ustream) for~x ∈ Γ (2.11)Combining potential flow with the Biot-Savart kernel seems at first glance tosatisfy requirements for adding boundary conditions to the free-space formulation.We have the Helmholtz-Hodge representation of a divergence-free vector field thatsatisfies our required boundary conditions. Unfortunately, the Biot-Savart kernelitself introduces errors into the velocity field. Following Saffman [64], if we re-compute the vorticity from a velocity field constructed by the Biot-Savart law, we10have∇×~u = ~ω+ 14pi∫Γ(~ω ·~n)∇ 1|~x−~x′|dS′ (2.12)The surface integral on the right of equation (2.12) is an error term that decaysby the inverse of the distance to any normal component of vorticity on the bound-ary. In order to correct this, one has to either introduce boundary conditions to thePoisson solve involved in creating the Biot-Savart law, or else employ a techniquewhere one adds fictitious vorticity in the interior of solid objects that creates a con-tinuous connection between the field lines of the normal vorticity. This fictitiousvorticity is known as image vorticity. In both cases, determining the boundary con-ditions or image vorticity is non-unique. For simple problems, however the lattermethod is probably the better technique.No Lagrangian vortex animation papers that have treated solid boundaries haveaddressed the problem of incident normal vorticity. Some methods, such as thosedescribed by Weissman et al. [74] and Brochu et al. [14], describe vorticity in termsof filaments or sheets. These don’t permit normal vorticity simply by requiring thesheets or filaments to be connected in the fluid volume. Normal vorticity wouldrequire the sheets or filaments to intersect the solid boundary.The neglect of incident normal vorticity isn’t completely restricted to the ani-mation community; Cottet and Koumoutsakos don’t mention it either in their book[19] on computational vortex methods.It could well be that for practical usage in animation, the visual nature of errorsgiven by not correctly treating incident normal vorticity are acceptable. In this case,the velocity field given by the erroneous formulation is still divergence-free andtherefore particle motion will still exhibit fluid-like motion. Also, on the molecularscale, fluids actually have a no-slip boundary condition; the no-through conditionis a subsequent approximation to the boundary layer motion of the fluid molecules.In this context, normal vorticity cannot exist, and the vorticity lines must passthrough the boundary layer to the far-side of the object in a connected fashion,which implies simulations like those in the papers by Weissman et al. [74] andBrochu et al. [14], limiting normal vorticity in the first place, are not too far off themark.112.3 Vortex Sheets and FilamentsOne of the crucial aspects of the vorticity equations is that they can be satisfiedby singular distributions of vorticity. These come in two types, vortex sheets andvortex filaments. Vortex sheets constitute a distribution of vorticity on a two di-mensional manifold. This distribution must still obey Helmholtz laws; the vortexlines must still be closed and incident on the manifold. As a result, for simplyconnected sheets, the vortex sheet is representable by a scalar potential, for whichthe vorticity follows its iso-contours.A common way to use vortex sheets is to place them on the boundaries of thesimulation and constrain them to cancel the normal flow. Park et al. [60] useda least squares method to constrain a vortex sheet on the boundary to cancel allincident flow on the surface. It’s not apparent that this technique maintained theunderlying characteristics of a vortex sheet or that it was guaranteed to enforce no-through boundaries. Weissman et al. [74] properly represented the sheet’s naturebut were required to perturb the resulting matrix to overcome the inherent null-space. Common to all these techniques is the fact that defining a vortex sheet tocancel no-through boundary conditions is simply equivalent to adding potentialflow into the domain. In chapter five of this thesis [14], potential flow boundarieswere implemented using an integral equation formulation that enabled a faster andmore efficient numerical method then the previous animation techniques we havejust discussed.12Chapter 3Boundary Integral EquationsBoundary integral equations were developed at the end of the 19th century. Theygave the first “proof” as to the existence of solutions to Laplace’s equation, andsome methods of calculating those. However, since these proofs and solutions de-pend on some functional analysis to be understood and aren’t in closed form, theyare less popular as a PDE introduction, and therefore occupy a somewhat obscurecorner of PDE theory in most mathematical departments. However, in terms ofnumerically calculating the solution to PDE’s, boundary integrals have enjoyed astronger position as they can reformulate a volumetric problem into one treatingonly the boundary or surface of the PDE domain. This dimensional reduction hasled to a variety of numerical treatments. Engineering literature often calls thesemethods “Boundary Element” methods.3.1 Boundaries and IntegralsThe “boundary in boundary integrals comes from the fact that we are solving aPDE by treating only the boundary of the domain. More specifically, we focus onLaplaces equation,∇2Φ(~x) = 0 ~x ∈Ω, (3.1)with a boundary condition specified for ~x ∈ Γ, the boundary of the domain Ω.For example, the interior Dirichlet problem prescribes the value of Φ on Γ with13Ω compact; the exterior Neumann problem prescribes the normal derivative of Φwith Ω the connected complement of a compact set (extending out to infinity).The basic assumption of our BIE theory is that we can construct a solution toLaplace’s equation by a distribution of charge on its boundaries and requiring thedistribution to satisfy boundary conditions. Suitably modifying equation (1.8) byrenaming our “charge” function f to σ and restricting it to a smooth boundary Γof a closed simple domain Ω givesΦ(~x) =∫Γ~yσ(~y)4pi|~y−~x|ds. (3.2)This type of distribution is called a single layer potential in integral equationliterature. A similar construction is called the double layer potential,Φ(~x) =∫Γ~yσ(~y)∂∂~n~y14pi|~y−~x|ds. (3.3)~n is the normal vector on Γ pointing away from Ω. It is easy to verify that both ofthese constructions solve Laplace’s equation in Ω. We simply apply the Laplacianto Φ, and bring it inside the integral. The trick is to now determine if all solutionson the domain can be represented in such a fashion. This requires determining themap of σ to the appropriate boundary conditions, which requires some tricky limittaking as we evaluate equations (3.2) and (3.3) closer to the boundaries. Both ofthe potentials are well defined for x ∈ Γ, but provide different benefits to the lim-iting process depending on which type of boundary condition is being solved for.Specifically, the single layer potential is used to solve the Neumann problem andthe double layer potential is used for solving the Dirichlet one. Our applications influid animation generally require solving the exterior Neumann problem.3.2 Integral EquationsFirst, we consider the interior Dirichlet boundary value problem, where we pro-scribe Φ(~x) = f (~x) on Γ and consider the solution on Ω. A first assumption is thatone could just evaluate Φ(~x) on Γ and arrive at the equation for σ by substituting14in f . For (3.2) this is exactly the case.f (~x) =∫Γ~yσ(~y)4pi|~y−~x|ds, for~x ∈ Γ. (3.4)However, for the double layer potential (3.3), we need to be more careful. Thefunction described by the double layer is actually discontinuous at the boundary Γ.Determining its value on the boundary requires careful limit taking and integration.The result isf (~x) =−σ(~x)2+∫Γ~yσ(~y)∂∂~n~y14pi|~y−~x|ds for~x ∈ Γ. (3.5)The extra term results integrating the discontinuity across the boundary. Thisterm is often referred to as a “jump condition”. In the case of integral equations,this jump condition is very desirable in proving the existence of solutions, in theconditioning of the problem, and in the ease of numerical treatment. We will touchon some of these points in the next section. First, though, we describe the exteriorNeumann boundary value problem. Here ∂∂~nΦ(~x) = g(~x) on Γ, and we remind thatthe connected domain Ω is defined by the complement of an appropriate boundedset, and thus includes infinity. Again, carefully take the single layer potential (3.2),and after similar limiting process, we arrive atg(~x) =σ(~x)2− ∂∂~n~x∫Γ~yσ(~y)4pi|~y−~x|ds ~x ∈ Γ. (3.6)Again, we have the jump condition, and the observation that the equation forthe exterior Neumann boundary condition differs from the integral equation for theDirichlet problem using the double layer potential simply by swapping ~x and ~y.(This means the two operators are adjoint.)153.3 Compact Operators and Fredholm TheoryWe can view our equations as linear operators acting on σ to give φ . These opera-tors are either of the formφ(~x) = I(σ) =∫ΓK(~x,~y)σ(~y))dy (3.7)orφ(~x) = (1+ I)(σ) = σ(~y)+∫ΓK(~x,~y)dy. (3.8)These two forms of equations are called Fredholm Integral equations of theFirst and Second Kind, respectively. Fredholm theory states that under some con-ditions on K, equations of the form (3.8) have either unique solutions, or a non-zero nullspace. While this statement is always true for equations involving finite-dimensional vector spaces, such as linear equations on Rn, it’s actually not neces-sarily true when dealing with infinite dimensional vector spaces, such as the spaceof square integrable functions on Γ. Equipped with this information, theorists atthe beginning of the twentieth century were able to prove existence and uniquenessfor Laplace’s equation by proving the same for σ in equations (3.5) and (3.6).Of critical importance to Fredholm theory is whether or not the operator∫Kσis a compact operator. We won’t give a rigorous definition of what a compactoperator is, as the functional analysis required is outside the scope of this thesis.However, we can offer some intuition and a definition for our context. Consider thekernel K on smooth bounded boundary Γ: it is compact if it is square integrable,meaning∫Γ∫Γ|K(~x,~y)|2d~xd~y < ∞ (3.9)and ∫Γ|K(~x,~y)|2d~x < ∞ (3.10)16and ∫Γ|K(~x,~y)|2d~y < ∞. (3.11)Determining the inverse for compact operators is provably ill-conditioned ifnot ill-posed. In the context of our integral equations, we are convolving over asurface. For a bounded smooth kernel, this convolution is a smoothing operation,and inverting a smoothing operator is in general ill-conditioned. Nevertheless, wesaw how convolving with the Dirac delta function in the introduction chapter is theidentity operation; obviously this convolution operator is not compact.Since one can define the Dirac function as the limit of a sequence of operatorswith bounded smooth kernels (see Griffel’s first chapter in his book on functionalanalysis [32]), we can conjecture that this ill-conditioning depends on the kernel.In fact, in Keeler [41], numerical experiments with a first kind equation using anintegrable singular kernel, the two dimensional analog of (3.4), showed only a lin-ear increase in condition number dependent on boundary discretization resolution.The critical part of these experiments was using high order interpolation whichefficiently captured the singularity.In comparison with first kind equations, second kind equations are much easierto treat numerically. The compact part of the equation can be regarded as a per-turbation of the identity operator, which allows for more regular computation, orregularization when the kernel is weakly singular. For example, Atkinson in hisbook, [6], proved a bounded condition number for a numerical solution to a secondkind boundary integral equation. In practice, when solving stable numerical lin-ear approximations to second kind equations, iterative methods often converge ina fixed number of iterations, independent of the problem size (albeit dependent onthe problem geometry). These characteristics of the second kind equations makethem far more desirable then the first kind, in both terms of efficiency and stability.3.4 Using Integral Equations for Potential FlowThe exterior Neumann problem (3.6), is uniquely suited to providing the poten-tial flow necessary for resolving no-through boundary conditions in vortex basedanimation techniques, as shown by Brochu et al. [14].17The equation has a unique solution even with multiply connected domains andimplicitly handles an infinite domain, which is also a strength of vortex methods.One drawback could be that it only provides a solution with zero circulation in mul-tiply connected domains. However, one can always define a circulation by loopingan image vortex filament around the boundary of a connection in a domain. Whilea vortex sheet can represent a circulation, each connection creates an additionaldegree of freedom in the representation that requires extra treatment for numericalsolutions.18Chapter 4The Fast Multipole MethodThe Fast Multipole Method (FMM), initially introduced by Greengard and Rokhlin[31], reduces the asymptotic complexity of the N-body problem from O(N2) toO(N) for particle-to-particle interactions arising from potential theory. We willgive a short description of the FMM with some details relevant to our implementa-tion, leaving most of the finer details to published literature.We consider the problem of computing the potential Φ and field ∇Φ given bya set of N point particles with charges qi and positions~xiΦ(~x) =N∑i=1qi4pi|~x−~xi| . (4.1)This leads to the N-body problem for finding the potential at each particle’s posi-tion,Φ(~xi) =N∑j 6=iq j4pi|~xi−~x j| . (4.2)Outside of the particles’ positions, the potential solves Laplace’s equation∇2Φ(~x) = 0, (4.3)Φ→ 0 as x→ ∞19The FMM uses the following series approximations to harmonic functions,in terms of spherical coordinates (r,θ ,φ) and the complex spherical harmonicsY mn (θ ,φ), to accelerate computing (4.2):Φ=∞∑n=0n∑m=−nFmnY mn (θ ,φ)rn+1, (4.4)Φ=∞∑n=0n∑m=−nLmnrnY mn (θ ,φ). (4.5)Far-field series expansions (4.4), also called multipole expansions, are realized bythe complex coefficients Fmn. These expansions satisfy that the far-field decaysat infinity and truncated series are used to efficiently approximate the potential ofclusters of particles evaluated at a distance from their location. Local series ex-pansions (4.5) with complex coefficients Lmn don’t satisfy the far-field conditionat infinity. However, the FMM uses truncated local expansions to accumulate ap-proximations to components ofΦ in a localized area due to distant particle clusters,akin to a Taylor Series. One can calculate the following far-field expansion for thepotential due to a group of point charges centered at the origin,Φ(~x)≈N∑i=0qi4pi|~x−~xi| =14pip∑n=0n∑m=−nFmnY mn (θ ,φ)rn+1(4.6)Fmn =N∑i=0qirni Ymn (θi,φi) (4.7)This approximation is valid for r > maxi |ri| and is truncated at n = p.The FMM commences by creating a hierarchical subdivision of particle posi-tions using an octree. The FMM then constructs the far field expansion for eachcollection of particles in each of the octree’s subdivision nodes containing a min-imum number of particles. This is done in O(N) time by computing the expan-sions from the bottom up and exploiting existing far-field expansions from childrennodes to construct that of their parent. We will refer to the calculation of far-fieldexpansions from existing ones as a Far-to-Far operation (FtF).Mild controversy surrounds the FMM’s claim to O(N) complexity, since straight-forward computation of the octree data structure required by the FMM is O(N logN).20In practice the creation of the tree is dominated by the approximation calculationsof the algorithm. It is actually straightforward to construct an octree in linear timeusing a binning or counting sort. In particular, we compute the octree in O(N)time using radix sort to create a linear octree. To do so, we specify beforehandthe maximum depth of the octree. We then compute unsigned integer coordinatesfor each particle that correspond to the node of the finest octree subdivision thatcontains the particle. A single key is then created by interleaving the bits of eachof the three coordinates, thus enumerating the nodes in Morton order [54]. We sortthe particles by these interleaved keys using a radix sort in O(N). The result of thesort is that the particles are arranged in a linear octree, from which the hierarchicaloctree structure can be constructed in O(N).Once the far-field expansions are computed, a downward pass is then con-ducted, in which local approximations at each octree node are used to accumulatethe potential due to adequately distant far-field expansions. This is again done inO(N) time by starting at the top level of the octree and proceeding downwards.Each node combines the approximation of the effect of relevant far-field expan-sions at the same level along with those local approximations already computed atthe parent node. The operation of computing local expansions approximating far-field expansions we refer to as Far-to-Local (FtL). Computing a local expansionfrom a parent node’s local expansion is called a Local-to-Local operation (LtL).The potential at each particle is finally calculated directly from neighbouringparticles that were never included in an adequately distant far-field approximationand a single computation of the finest possible local approximation. The fast-multipole method as stated is provably O(N) for an even distribution of particles;adaptive techniques also exist for uneven distributions. The asymptotic optimalityof the FMM comes at a cost of significant overhead which is O(p4)with the trunca-tion order p due to the matrix representation of the translation operations FtF, FtL,and LtL. It was discovered that these could be reduced to O(p3) by rotating coeffi-cients into a desired coordinate frame before computing the operations, which ourimplementation in chapters 5 and 6 utilize. This high overhead cost has dampenedenthusiasm for the three dimensional FMM in in the scientific computing commu-nity. However, due to the relaxed need for accuracy in fluid animation, we usep = 4 as a standard, and the FMM quickly overcomes its asymptotic constant and21becomes relevant for representative particle counts, as we will see in the results ofour next section.22Chapter 5Linear-time Smoke Animationwith Vortex Sheet MeshesThis chapter contains the paper “Linear-time smoke animation with vortex sheetmeshes” by Brochu, Keeler and Bridson [14]. The explicit mesh tracking and vor-tex sheet implementation were contributed by Brochu while the Fast MultipoleMethod implementation and potential flow boundaries formulation based on inte-gral equations were Keeler’s. The boundary conditions in particular solved a seri-ous problem with pre-existing Lagrangian vortex methods published in animationliterature.5.1 SynopsisWe present the first quality physics-based smoke animation method which runs intime approximately linear in the size of the rendered two-dimensional visual detail.Our fundamental representation is a closed triangle mesh surface dividing spacebetween clear air and a uniformly smoky region, on which we compute vortexsheet dynamics to accurately solve inviscid buoyant flow. We handle arbitrarymoving no-stick solid boundaries and by default handle an infinite domain. Thesimulation itself runs in time linear to the number of triangles thanks to the useof a well-conditioned integral equation treatment together with a Fast MultipoleMethod for linear-time summations, providing excellent performance. Basic zero-23albedo smoke rendering, with embedded solids, is easy to implement for interactiverates, and the mesh output can also serve as an extremely compact and detailedinput to more sophisticated volume rendering.5.2 IntroductionTantum videt, quantum computato: as much as one sees, that much one shouldcompute. While the last decade and a half has seen a revolution in the quality ofcinematic smoke animation through the use of physical simulation, scalability re-mains a serious impediment to high quality real-time smoke effects and interactiveartistic design. To achieve ∼ n× n visual detail in 2D renders, all prior methodshave needed at least O(n3) time per step, due to the simulation or the rendering orboth. In this paper we develop a purely Lagrangian vortex sheet model of smokesimulation, and argue that in the most useful asymptotic limit it captures all essen-tial dynamic and visual detail in a closed 2D surface. We discretize the model withan adaptive dynamic triangle mesh, implement time integration using an optimallinear-time Fast Multipole Method (FMM), and render the result at interactive rates— computing only what is seen. Additionally, our method handles arbitrary no-stick moving solid boundaries and provides output for more sophisticated offlinerendering.Previous smoke simulation methods in graphics can loosely be divided betweenvelocity-pressure formulations and vorticity approaches. The classic example of avelocity-pressure solver is Stam’s Stable Fluids approach [69]: this requires an n3grid of the volume, and even with an optimal linear-time multigrid solver for thepressure projection, is an order of magnitude more costly than the desired O(n2)rendered output. Vortex methods use a much richer representation whereby inprinciple many fewer elements (such as our vortex sheet of O(n2) triangles) canreplicate the same detailed dynamics. However, so far in graphics the solvers hit asimulation bottleneck when calculating velocity from the vorticity and solid bound-ary conditions, resulting in a dense n2× n2 matrix solve which can be as bad asO(n4) storage and O(n6) run-time using LU factorization. Prior vortex methodsalso have used simpler volumetric soot particle tracking for providing output tothe renderer, which introduces another expensive O(n3) bottleneck. For film visual24Figure 5.1: We represent smoke with an adaptive triangle mesh both forlinear-time simulation, as a vortex sheet, and linear-time interactive ren-dering as the boundary of the smoky region. Here a deforming torusinfluences the buoyant smoke flowing around and through it.effects, this strains the file system and network capacity in particular, as simulationis generally run separately from rendering, and so all particle data must be storedand transferred to the renderer.Our basic model in this paper treats air as an incompressible inviscid fluid withthe Boussinesq buoyancy approximation,25∂~u∂ t+~u ·∇~u+ 1ρ∇p = α~g (5.1)∇ ·~u = 0 (5.2)where α~g is the buoyancy force, with α = 0 in the ambient clear air, and α in thesmoky region depending on the ratio of smoke temperature to ambient temperatureas well as the density of soot particles. At solid boundaries we take the no-stickboundary condition,~u · nˆ =~usolid · nˆ, (5.3)and at infinity the far-field condition ~u = 0. We assume ~u = 0 initially. Withoutviscosity, the vorticity ~ω = ∇×~u remains zero apart from where buoyancy createsit; vorticity is advected and stretched by the flow, but doesn’t diffuse away fromsmoke/temperature gradients.Just as the viscosity in air is virtually negligible for most practical scenarios,hence the inviscid equations are an excellent model, the rate of molecular diffusionof temperature and soot particles in air is vanishingly small, so they are just ad-vected through the velocity field. Under the assumption that smoke sources haveuniform temperature and soot concentration, and also using the incompressibilitycondition, the air is precisely divided into uniformly clear and uniformly smokyregions. This in turn implies buoyant vorticity is only generated at the surfaceseparating the two, where α jumps in the normal direction, hence all vorticity isconcentrated in that surface, making a vortex sheet. As Stock et al. show [70], thevorticity is further constrained to remain tangent to this sheet — in section 5.4, weshow that our discretization identically enforces this constraint.5.3 Related WorkFluid simulation for computer animation has become a broad topic in recent years.A thorough survey of techniques is beyond the scope of this paper, so we referthe reader to the textbook by Bridson [9], and focus on work related to smoke.Foster and Metaxas animated smoke using a one-phase, grid-based Navier-Stokes26solver [26]. They used marker particles to track smoke through the simulation,transferring particle density onto a grid at render time to make use of volumetricrendering techniques. Stam introduced an unconditionally stable time integrationscheme to computer graphics [69], and his method was later adapted to use a stag-gered grid with vorticity confinement [23]. His method uses a grid-based scalardensity field to track the concentration of smoke. Recent work has focused on re-ducing the damping effect of numerical diffusion by introducing improved integra-tion schemes [44, 67, 83], adding sub-grid-scale turbulence on top of a simulation[46, 57, 65], or through the use of procedural methods [10, 66].Losasso et al. improved the scalability somewhat with the replacement of theuniform grid by an octree structure, but adaptivity in velocity-pressure formulationis difficult as much of the volume must still be finely gridded for adequate be-haviour [52]. Other fluid solvers (e.g. [22]) only run a solve on a sparse tiled gridinstantiated in a bounded envelope around the smoke, but need a substantial enve-lope to approximate an unbounded domain, and run at full volumetric resolutioninside the smoke, likewise limiting scalability.The use of triangulated surface meshes in fluid animation has traditionally beeneschewed in favour of grid- or particle-based surface tracking methods, such asthe particle level set method [21]. Recently, however, some researchers in thecomputer graphics community have begun turning to triangle mesh-based surfacetracking for liquid simulation [13, 55, 75, 76, 80]. These authors generally useexplicit surfaces to track the fluid interface in a free-surface simulation, whereas inthis paper, we use surfaces themselves as simulation elements.Mesh surfaces have been used for visualizing 3D vector fields since Hultquistgeneralized streamlines to stream surfaces [36]. More recently, streak surfaces(moving open surfaces seeded continuously from a curve), and time surfaces (seededinstantaneously from a surface) have been treated using surface tracking methodssimilar to the one used in this paper. Krishnan et al. [47] recently introduced analgorithm which advects an adaptive mesh through a 3D vector field to treat bothof these categories of surface. A similar visualization method which explicitly usesa smoke metaphor was introduced by von Funck et al. [73]. Their method does notuse adaptivity, but instead gradually reduces the opacity of poor quality triangles,giving the appearance of dissipating smoke or steam.27Brochu & Bridson [12] used surface tracking to visualize a procedural smokesimulation. Park et al. [59] constructed NURBS surfaces from particle streamlinesto achieve thin, sheet-like surface renderings of smoke simulations. These tech-niques are useful for visualizing an external simulation, however neither approachuses the surface elements for simulation, as we propose in this paper.Vortex methods for fluid dynamics were developed during the 1980s as an ef-ficient and high-order, purely Lagrangian method of numerically solving the Eulerequations for fluid flow [3, 8]. Subsequent work was done on developing vis-cous vortex methods for the Navier-Stokes equations and semi-Lagrangian con-texts [19]. Though the grid distortion in fully Lagrangian vortex methods makesthem less useful for scientific application, their efficiency in preserving rotationalflow makes them appealing for fluid animation.In computer graphics, early 2D simulations discretized vorticity on a grid [79],or used vortex-in-cell methods, advecting Lagrangian parcels of vorticity throughthe domain but relying on a grid-based Eulerian solver [18, 28]. Park and Kim[60] introduced a fully Langrangian vortex particle method, using the Biot-Savartformula to compute velocity on the particles without a grid. Recently, Weißmannand Pinkall [74] have demonstrated a complete vortex filament method, with adap-tive filaments. Both approaches discretize solid boundaries with vortex sheets andrely on static geometry to enable precomputation for efficiently handling boundaryconditions. In comparison, our linear-time approach can handle arbitrary movingsolid boundaries, as demonstrated in section 5.5.Several papers have used vortex primitives to augment Eulerian fluid simula-tions, including vortex particles [66] and Eulerian vortex sheets [45]. Concurrentwith this work, Pfaff et al. have introduced a Lagrangian vortex sheet method foraugmenting an Eulerian fluid solver [62]. Their vortex sheet discretization is sim-ilar to the one presented in this paper; however their vortex sheet interactions arepurely local, with global flow handled by an underlying grid-based, Eulerian fluidsolver. Their use of compact kernels avoids the need for FMM, but the reliance onthe grid-based solver limits scalability and complicates the simulation process.Boundary integral equations were initially developed a century ago as a tool forproving the uniqueness and existence for solutions to PDEs, particularly Laplace’sequation. They have had a modern renaissance as effective numerical methods for28solving a variety of PDEs [30]. This has been primarily motivated by the con-struction of effective summation methods such as the FMM which allows rapidsolution of the dense linear equations resulting from the numerical treatment ofmany integral equations.5.4 The Vortex Sheet MeshOur discretization of the vortex sheet closely follows Stock et al. [70] and Pfaff etal. [62] (though unlike their methods which rely on a 3D background grid for bulkvelocity computation, we directly compute it in linear time on the surface itselfwith FMM, including arbitrary solid boundaries, as we will see in section 6.4). Inparticular, we discretize the vortex sheet surface with a triangle mesh, and storecirculation values, Γi, on each mesh edge i. The per-triangle vortex sheet strengthis then:~γp =1ap∑j~e jΓ j, (5.4)for triangle edges j, where ~e j is the vector between the end points of edge j, andwhere ap is the area of triangle p.The vorticity on a triangle is then computed as required with the formula:~ωp = ap~γp =∑j~e jΓ j. (5.5)for triangle edges j, where ~e j is the vector between the end points of edge j.The main advantage of using these scalar circulations as our fundamental variableis that they are conserved along particle paths. The vorticity stretching is thusimplicitly handled by the stretching of the vortex sheet mesh. Recomputing thevortex sheet strength when needed for the velocity update is done according to(5.5). Using (5.5) also guarantees that ~ω is tangent to the triangle plane.We apply a buoyancy force α~g = 〈0,α,0〉 to each triangle as ∆~ωi = α~g×~niAifor triangle normal ~ni and area Ai. We then compute the change in circulation foreach triangle edge by solving the overdetermined system of equations 5.5 for Γ assuggested by Stock et al. [70]. We solve this in the least-squares sense. We notethat recomputing the total edge circulations in this way would introduce numeri-29cal damping manifesting as undesirable loss of vorticity, so we only solve for thechange in edge vorticity due to the buoyancy force.With no solid boundaries, the velocity at any point in the domain can be com-puted explicitly with the Biot-Savart law, with S being the vortex sheet surface, andK being a (mollified) fundamental solution to Laplace’s equation:~ufluid(~x) =∫S∇K(~x−~y)×~ω(~y)dy (5.6)This is discretized on a mesh surface using midpoint quadrature as:~ufluid(~x)≈∑j∇K(~x−~c j)× ~ω j (5.7)where~c j is the barycentre of triangle j, and ~ω j is the constant vorticity on trianglej.Since we require the velocity at each vertex, and every velocity evaluation de-pends on the entire surface, this direct summation has computational complexityO(M2) in the number of mesh vertices. Later in this paper we will discuss usingthe Fast Multipole Method to achieve an approximation which runs in O(M) time.The gradient of the fundamental solution ∇K(~x−~y) is singular when~x=~y andbecomes large as~x approaches~y, therefore a mollified approximation is often used.Our implementation follows Beale and Majda’s form [7]:∇K˜(~x−~y) =−(~x−~y)1− e(r/β )3r3, (5.8)where r = ||~x−~y||, and β is a mollification parameter.As our mesh evolves due to this self-imparted velocity, it naturally expands andcurls up. Without some form of surface mesh refinement, the simulation wouldquickly lose accuracy and visual appeal. On the other hand, we wish to keep thenumber of mesh elements as low as possible, to keep the computational cost down.Therefore, removing small elements (recoarsening) is also desirable. For meshrefinement and recoarsening operations, we use the open-source El Topo surfacetracking library [11]. This library also has the ability to perform changes in topol-ogy, such as merging and pinching. By merging nearby surface patches, we can30eliminate very thin sheets, thereby further reducing the total number of mesh trian-gles while still maintaining a good quality simulation.5.4.1 Circulation RedistributionEl Topo’s surface remeshing operations add and remove edges. As we are storingscalar circulation values on the edges, we must update these values as the meshconnectivity changes. One way of doing so is to compute per-triangle vorticity val-ues using equation 5.5, performing the remeshing operations, then solving equation5.5 for Γ in a least-squares sense. However, we have found that repeatedly solvingthis overdetermined system leads to diffusion of vorticity.An alternative is to explicitly update the edge circulation values using somerules of thumb when a remeshing operation is performed. Though heuristic, wehave found these rules reduce vorticity diffusion as compared to repeatedly pro-jecting back and forth between edge and triangle values. Of the local remeshingoperations available in El Topo, we use only edge splitting, edge collapsing, andmesh merging via deletion and creation of new triangles.Edge splitAs shown in figure 5.2, splitting edge AB with opposite vertices C and D introducesa new vertex E. We store the circulation value on the original edge AB, and assignit to both of the new edges:ΓAE = ΓEB = ΓAB. (5.9)The two other new edges, CE and ED, receive no circulation. This ensures thatthe sum of vorticities over the new triangles (by equation 5.5) equals the sum ofvorticities on the new triangles:~ωABC+~ωDBA = ~ωAEC+~ωCEB+~ωDBE+~ωDEA (5.10)31Figure 5.2: Edge split operation. New edges AE and EB are assigned circu-lation from original edge AB.Edge collapseWhen collapsing an edge, we simply allow the circulation on that edge to disappearfrom the system. Since we collapse only very short edges, the loss in circulation isgenerally minimal.Mesh mergeEl Topo merges surface meshes by removing two nearby edges and their incidenttriangles, and introducing new triangles, forming a neck between the surfaces [11].From the new edges introduced, we pick the two edges whose associated vectors(difference of end points) best match the removed edges, i.e. finding new edges iand j such that||~ei−~e0||+ ||~e j−~e1||is minimized, where edges 0 and 1 are the removed edges. We assign the originaledge circulations to these “best match” new edges (see figure 5.3).5.5 Solid BoundariesBoundary conditions for Lagrangian vorticity-based fluid animation have generallybeen constructed by defining a vortex sheet on the surface of the solid boundaries.Park and Kim [60] implemented boundary conditions in 3D by enforcing the tan-gential and normal components of the fluid velocity to be equal to those of the32Figure 5.3: Edge merge operation. Left: The red and green edges and theirincident triangles are removed from the mesh, leaving two quad-shapedholes. Right: The resulting holes are zippered with new triangles. Sup-pose the red and green edges in the right diagram are the closest matchesto the corresponding red and green edges in the left diagram (comparingEuclidean distances of associated edge vectors). These new edges arethen assigned the circulation values from the deleted edges.boundary using a vortex sheet, which resulted in an over-determined matrix solve.To avoid finding an optimal solution each time-step, they restricted themselvesto a simple static boundary for which they precomputed a pseudo-inverse, requir-ing only a matrix-vector multiplication instead of a matrix solve. Weißmann andPinkall [74] defined their vortex sheet in terms of a scalar potential, avoiding theextra dimensionality of the tangential vorticity field. They regularized the resultingmatrix and also precomputed a factorization, again limiting their domain to a singleboundary with static geometry.We incorporate solid boundaries by adding an irrotational, divergence-free vec-tor field that corrects the boundary velocity induced by the vorticity to account forsolids. To do so we define the irrotational flow by the gradient of a scalar potentialfunction Φ that satisfies Laplace’s equation ∇2Φ= 0 in the domain. The potentialdue to a single point charge with value σ , at point~y, is:σ4pi|~x−~y| . (5.11)We represent Φ by a continuous distribution of point charges {σ(~y)|y ∈ S}33where S is the boundary of our domain.Φ(~x) =∫Sσ(~y)4pi|~x−~y|dy. (5.12)This representation of Φ is called the single layer potential. The normal derivativeof the single layer potential has a discontinuity across the domain boundary. Whensatisfying boundary conditions, this leads to the following Fredholm integral equa-tion of the second kind for σ [35, 43]. Given the normal derivative f of Φ on theboundary we have:f =−σ(~x)/2+∫Sσ(~y)∂∂nx14pi|~x−~y|dy, ~x,~y ∈ S. (5.13)Fredholm equations of the second kind are integral equations consisting of theidentity operator plus an integral operator that is “compact”. While the definitionfor compact requires an understanding of functional analysis and is therefore be-yond the scope of this paper, an integral operator is compact if the multiplicativeterms in the integral are themselves square integrable, which is the case in (5.13).The solvability of Fredholm equations of the second kind can be determined byexamining the null-space of the equation, similar to determining the solvability oflinear systems. For smooth boundaries, (5.13) has only a trivial null-space andthus its solution exists and is unique for any right hand side f . In addition, solvableFredholm equations of the second kind tend to lead to well conditioned numer-ical approximations. While our description of integral equations is cursory andincomplete, there exists ample mathematical literature for the interested. Solving(5.13) for σ and evaluating Φ from equation (5.12) satisfies the exterior Neumannproblem for Laplace’s equation,∇2Φ = 0 (5.14)∂∂nΦ(x) = f , x ∈ S. (5.15)Let ~ufluid be the velocity due to the fluid, and ~usolid be the velocity of the solidobject. Setting f = (~usolid−~ufluid) ·~n, solving for potential Φ, and calculating thegradient yields a potential flow which will cancel the normal component of velocity34at the solid boundary, effecting a free-slip, no-entry boundary condition.The advantages of this scalar formulation over the previously described vortexsheet boundary implementations are that it is solvable for arbitrary smooth solidboundaries and is particularly amenable to numerical solutions based on iterativematrix solvers and fast summation methods. Integral equations for implementingboundaries based on vortex sheets are more difficult to implement, and though onecan derive Fredholm integral equations of the second kind for finding a bound-ary vorticity [19], these formulations become singular if the domain is multiply-connected.Discretizing (5.13) with midpoint quadrature yields the following equation fortriangle i:fi ≈−σi/2+∑jσ j∂∂ni14pi|~ci−~c j|A j (5.16)Directly evaluating this summation for all triangles yields an M× M linearsystem where M is the number of triangles on the solid boundary. In practice,this matrix is very well-conditioned; for such problems, iterative Krylov solverssuch as BiCGSTAB have been observed to converge in O(1) iterations, resulting intotal complexity of O(M2) for the solid solve. We can additionally accelerate thisto O(M) by using the FMM to compute the matrix-vector product in the iterativesolve rather than explicitly constructing the system.Once we have solved for a density σ , we can compute Φ as in equation (5.12).The gradient of Φ is then evaluated as:∇xΦ(~x) =∫S−σ(~y)4pi|~x−~y|3 (~x−~y)dy (5.17)The midpoint quadrature rule for this integral yields:∇xΦ(~x)≈∑j−σ jA j4pi|~x−~c j|3 (~x−~c j) (5.18)We use this to modify our fluid surface vertex velocities:~ufinal(~x) =~ufluid(~x)+∇Φ(~x) (5.19)355.6 Fast Multipole MethodThe Fast Multipole Method (FMM), initially introduced by Greengard and Rhoklin[31], reduces the asymptotic complexity of the N-body problem from O(N2) toO(N). The FMM is used in two different areas of our fluid simulation. Computingthe velocity from the vorticity using the Biot-Savart law can be accomplished usingan FMM for each component of the vorticity. In addition, computing the matrixmultiply in the iterative solver for the solid-boundary potential flow can be reducedfrom O(M2) to O(M) by using the FMM to approximate Φ:Φ(~xi) =∑jσ j|~xi−~x j| , i, j = 1..N, (5.20)and its gradient,∇Φ(~xi) =∑j−σ j(~xi−~x j)|~xi−~x j|3 (5.21)for a set of points ~xi and “charges” σi. Our 3d FMM implementation follows thatdescribed by Cheng et al. [15].Mild controversy surrounds the FMM’s claim to O(N) complexity, since straight-forward computation of the octree data structure required by the FMM is O(N logN).We compute our octree in O(N) complexity. To do so, we specify beforehand themaximum depth of the octree. We then compute unsigned integer coordinates foreach particle that correspond to the grid cell of the finest octree subdivision thatcontains the particle. A single key is then created by interleaving the bits of eachof the three coordinates creating a 3D Morton ordering of the occupied grid cells.We sort the particles by these interleaved keys using a radix sort in O(N). Theresult of the sort is that the particles are arranged in a linear octree, from which thehierarchical octree structure can be constructed in O(N).The log-log plot in figure 5.4 shows the order of complexity for computingvelocity due to vorticity using the direct summation method vs. the FMM. Perfor-mance plots for the solid solver are shown in figure 5.5.36Figure 5.4: Simulation time per frame for computing velocity via the Biot-Savart law, for both the direct method and Fast Multipole Method.5.7 RenderingOur interactive-rate smoke rendering is accomplished through a simple set of OpenGLshaders, which compute the optical depth of the smoke volume. We then apply azero-albedo absorption model based on this depth, taking into account the solid ob-jects in the scene. More sophisticated real-time smoke rendering methods do exist(cf. [82]), and exploring self-shadowing and light scatter for triangle mesh modelsin real time is an area of potential future work.As Brochu and Bridson [11] point out, for high quality cinematic smoke ren-dering in a typical film pipeline, the critical bottleneck is often the load on the filesystem and network from tracking smoke per frame with hundreds of millions ofsoot particles, or similarly high resolution grids. The surface mesh representation isvastly more efficient for capturing the volumetric bulk of the smoke, and as Brochuand Bridson note can be directly ray-traced for self-shadowing smoke rendering.37Figure 5.5: Simulation time per frame for computing the single layer poten-tial on a sphere using BiCGSTAB, computing the matrix multiply withboth the direct method and Fast Multipole Method. Note that the di-rect method is nearly exactly O(N2) showing the O(1) convergence ofBiCGSTABWe extend their idea noting that the interior of the surface mesh, at render time, caneasily be rasterized into a lower resolution 3D texture for scattering computation,at least using a box filter via polygon clipping; if this grid is aligned with the lightsource or the view, especially efficient calculation of single-scattering is possible.The direct render can use triangle rasterization of the mesh, storing depths in anA-buffer, looking up into the 3D texture for smooth lighting calculations. Goingeven further, we can track the age of mesh elements since their emission from asource, and splat this into another render channel to provide age-proportional blur-ring, thus giving the visual effect of smoke dissipation from turbulent mixing ordiffusion at render time, for every frame in parallel. Figure 5.6 illustrates the effectof an image-space age-dependent blur.38Figure 5.6: Our proof-of-concept volume renderer produces, from left toright on the top row, an alpha channel from the smoke mesh, a scatteringchannel with self-shadowing, and an image-space age channel. Belowon the left is a sharp composite of the first two, and below on the rightincludes an age-dependent blur, simulating diffusion in image-space atrender-time.395.8 ResultsOur FMM-accelerated vortex sheet solver was implemented in C++, using El Topofor surface tracking, and tested on several Linux and Mac OS X computers. Figures5.4 and 5.5 show that our tests have empirically near-linear time complexity. Sincewe do not rely on precomputing the solution to the exterior Neumann problem, wecan easily handle non-rigidly deforming solid objects, as shown in figure 5.1.Our interactive smoke renderer implemented in GLSL runs at 4 FPS for 168Ktriangles on an AMD Radeon 6770M with 512 MB of memory (shown in figure5.1). Our proof-of-concept software renderer uses CPU-based volumetric rasteri-zation and self-shadowing to produce the images shown in figure 5.6. These imageswere produced from 126K triangles at 1600×1200 resolution with a 167×125×82 shadow texture in around 5.5 s on a single thread of a 2.3GHz Intel Core i7CPU, including all file IO.5.8.1 Mesh SimplificationAs mentioned in section 5.4, remeshing operations are crucial for allowing the sur-face mesh to move freely while keeping the number of triangles manageable. Theeffectiveness of edge splitting and collapsing to control the explosion of mesh el-ements has been shown by Stock et al. [70] and others, but we also make use ofEl Topo’s topology change operations — specifically the merging and pinching ofmesh surfaces. We have found this to be an effective tool in mitigating the explo-sion of surface area (and number of mesh elements), which is a known problem invortex sheet simulations.As an example, we ran a simulation twice, once with topology changes en-abled, and once with these changes disabled. The number of triangles per timestep is shown in figure 5.7. Note that without topology changes enabled, the num-ber of triangles grows exponentially, but with aggressive topological merging andsplitting, the number of triangles remains bounded (see the accompanying video).Pfaff et al. [62] also address this explosion in the number of mesh elements.Their approach identifies triangles which are deep within the volume of smoke(using a grid-based signed distance field), or otherwise occluded from the cameraview, and marks them for deletion. By contrast, our approach uses only Lagrangian40Figure 5.7: Geometry creation when simulating a smoke plume withouttopological changes (Without TC) and with topological changes (WithTC).mesh-based operations, without additional structures or heuristics.5.9 Future WorkWe have demonstrated that interactive-rate, near-cinematic-quality smoke simula-tion and rendering is within reach, but there are many avenues to explore further:• Code optimization and parallelization of the FMM solver to achieve interac-tive simulation rates.• Porting our software A-buffer rasterizer to run self shadowing and diffusionin GPU hardware.• View-dependent, level-of-detail mesh refinement, coarsening, and topologychanges.• Incorporating the creation of vorticity from solid interactions (vortex shed-ding).• Simulation of flame front propagation to achieve vortex sheet based fire sim-ulation, combined with real-time lighting techniques on the rendering side.41• Simulation of ocean wave free surfaces — much of our vortex sheet workcould apply to rough and vast simulations, which are currently infeasiblewith volumetric simulation methods.AcknowledgementsThis work was supported in part by a grant from the Natural Sciences and Engi-neering Research Council of Canada. We also thank Mary Bridson for the Latinturn of phrase.42Chapter 6Ocean Waves Animation usingBoundary Integral Equations andExplicit Mesh TrackingThis chapter contains the paper “Ocean Waves Animation using Boundary Inte-gral Equations and Explicit Mesh Tracking” by Keeler and Bridson [42]. This pa-per provides the first animation technique for deep ocean waves based on integralequations.Figure 6.1: Deep Ocean Waves6.1 SynopsisOcean waves are a common animation challenge. To achieve natural-looking andrich detail, physical simulation of one kind or another has generally been adopted.43However, there is an inconvenient gap between highly efficient but limited FFT-based surface displacement methods and full 3D simulation of a large volume ofwater, which we aim to bridge in this paper.For non-overturning waves without significant interaction with boats or othersolids, or just to provide a good ripple-scale animated displacement texture,Tessendorf’s extremely efficient FFT-based solver [71] is typically used. Thismodel is derived from an irrotational flow assumption (i.e. that the flow beneaththe surface has negligible vorticity) which has been experimentally validated forgentle to moderate waves, along with several other simplifications. However, itsrealism falters for stormy oceans, where the assumption of a small and smooth per-turbation from flat geometry is grossly broken, and doesn’t account for the presenceof large solids in the water.For full-fledged interaction with solids and for overturning waves with splashes,3D simulation of the free-surface Navier-Stokes equations is the main alternative.However, this comes with a significant performance cost, as a large volume of wa-ter beneath the surface must be discretized and solved: though it is not directlyrendered, the effect of the depth is plainly visible in terms of wave speeds and dis-persion. The finite simulated domain generally has to be much smaller than thefield of view in rendering as well, requiring nontrivial effort to convincingly blendthe 3D region into some sort of continued surface geometry out to the edges of theview.Our research provides a novel 3rd mode of wave animation based on rephras-ing the potential flow problem as a boundary integral equation. This allows allkinematic and dynamic properties of the waves to be represented and solved on thefree and solid boundaries of the fluid itself.The core ideas of our research parallels those stated in papers by Xue et al.[78]and Liu et al.[50] which present a boundary-integral formulation for deep oceanwaves and solid boundaries. Additionally, there also exists an FMM acceleratedwaves simulation [25]. Our contributions lie in the combination of these ideaswith adaptive and robust explicit mesh tracking to provide a novel tool for waveanimation. We also use our own unique formulation that is more amenable to usingthe FMM and allows analytic calculation of the flow gradients necessary to movethe surface without finicky surface geometry approximations.44Our method differs from those referenced in that we use an indirect integralformulation which obviates the need for high-order geometry and function approx-imation (required by the engineering simulations). As a result we are able to em-ploy explicit mesh tracking using the El Topo library [11] to represent the mesh ina robust and adaptive manner.There also exists an accelerated BEM wave simulation, developed by Fochesatoet al. [25], using principles based on the Fast Multipole Method (FMM), an accel-eration algorithm due to Greengard et al. [31]. Our reformulation is immediatelyamenable to the FMM and we use it to demonstrate the favourable asymptoticcomplexity of our simulation.In short, our method presents the following novel contributions to the anima-tion of water wave phenomena:• a surface-only memory-efficient method for computing deep ocean waves,• a boundary integral formulation better fitted to coarse explicit meshes, in-cluding solid boundaries,• a method which implicitly satisfies deep water conditions without needing alarge simulation domain,• optimal scaling of computation time using advanced hardware and algo-rithms.6.1.1 Related WorkThe use of BIE methods in graphics began with James and Pai [38] for elastic de-formation. However, the most common use has been in the computation of solidboundaries for vortex methods. Many methods construct sheets of singular vor-ticity on solid surfaces which then cancel normal flow. Solving for the necessaryvortex sheet strength constitutes solving a BIE: see Park and Kim [60], Weiss-man and Pinkall [74] and Golas et al. [29] for this approach to incorporating solidboundaries in vortex-based smoke simulations. Brochu et al. [14] used a differ-ent BIE formulation called the single layer potential to compute the potential flowsatisfying free-slip boundary conditions necessary in a vortex method, and demon-45strated using the FMM combined with iterative methods to accelerate computingthe solution.Wave animation in general has a long history in graphics, and we do not attemptto document it all. Looking at surface models which can be derived from potentialflow, Fournier and Reeves’ [27] and Peachey’s [61] physically-based proceduralmodels are a major milestone, using Gerstner or Rankine waves and also includingthe effect of solid boundaries such as beaches. Tessendorf’s FFT simulation [71]mentioned above also belongs to the potential flow / surface category. Other re-searchers have turned to discretizing the shallow water equations, which allows fora very efficient surface-based simulation [4, 16, 48] but only approximates shallowwaves: the realistic dispersion of deep water ocean waves is not captured. Yuksel’swave particle approach [81] similarly does not capture the dispersion of deep waterwaves.Full three-dimensional water solvers can be applied to ocean scenes, thoughusing a volumetric grid deep enough to effectively model deep water waves, broadenough to capture a reasonably large surface area, and detailed enough to get ad-equate surface resolution can make this extremely expensive. Irving et al. [37]and Chentanez and Mu¨ller [17] incorporated vertically stretched grid cells in theirsolvers to make deep water more economical. Nielsen and Bridson [58] guided ahigh resolution simulation of just the surface layer of the ocean with a coarse butmuch deeper simulation to capture deep water effects.6.2 Deep Water WavesThe key assumption in our non-linear waves formulation is that the wave motionis irrotational and thus can be modelled by potential flow. This assumption iscommonly held, and has experimental evidence supporting it. Potential flow is adivergence-free velocity field~u obtained by taking the gradient of a scalar potentialfunction Φ(~x), where Φ(~x) solves Laplace’s equation,∇2Φ= 0. (6.1)46Potential flow is uniquely defined by its boundary and the conditions thereon. It isnatural then to restrict our mathematical and computational model to the boundaryusing integral equations.Figure 6.2: The geometry of our wave formulation.For the deep ocean, as seen in Figure 6.2, our scalar potential exists in an in-finite domain Ω with boundary Γ comprised of solid boundaries Γs and the freesurface representing the water-air interface Γ f . We will denote points on the sur-faces with subscripts: ~x f ∈ Γ f ,~xs ∈ Γs and~xΓ ∈ Γ.On the free surface the solution to Bernoulli’s equations (which govern poten-tial flow) gives a time dependent Dirichlet condition,DDtΦ(~x f ) = 12 |∇Φ(~x f )|2−g z(~x f )− p(~x f ), (6.2)for the atmospheric pressure p, height z, gravitational acceleration g, and the con-vective derivative D/Dt = ∂/∂ t +(~u ·∇). At solid boundaries Γs in the fluid weimpose the normal velocity~n ·~us of the solid as a Neumann boundary condition,∂∂nΦ(~xs) =~n(~xs) ·~us(~xs). (6.3)To approximate the near infinite depth of the ocean floor, we apply the far-fieldcondition Φ→ 0 as z→−∞.The wave motion is given by advecting the surface Γ by the gradient of thepotential,47DDt~xΓ = ∇Φ(~xΓ) (6.4)The entire system is determined at a time t by the variables Φ(~x f ), Γ(t) andsolid boundary velocity. The system is initialized by defining these variables attime t = 0. We generally operate under the assumption that the atmospheric pres-sure p is zero on the boundary. Added modelling could include variable pressuresuch as that due to wind forces.6.2.1 LinearizationThe non-linear term in Bernouli’s equation (6.2) is measurably small for oceanwaves. We moreover find it can be destabilizing if included in the numerics,thus we linearize the equations by dropping the non-linear term, consistent withTessendorf’s commonly-used FFT wave simulations. We are similarly able to sim-ulate the shallow troughs and steep peaks of ocean waves without it, as we followthe full motion of the surface and not just the vertical displacement. We expect infuture work that more advanced numerical treatment will allow stable inclusion ofthis term when it is relevant, allowing for a wider range of animated phenomenawhere the non-linear term is more important such as wave breaking and overturningwaves.6.2.2 Boundary Integral FormulationThe potential flow wave formulation above lends itself directly to reformulation us-ing a boundary integral equation. Boundary integrals represent the solution interiorby integrating over data contained on its boundaries; finding the solution requiresmatching the boundary data to the given boundary conditions. Integral equationscome in two types, direct and indirect. In direct integral equations the boundarydata is some physically relevant quantity, such as the solution values or normalderivative on the boundary. The aforementioned engineering literature adopts thisapproach in their approach to wave modeling, using Green’s second identity torepresent the interior solution as an integral involving both Dirichlet and Neumanndata. For reasons below, we instead take an alternative approach, an indirect inte-48gral formulation, where boundary data used to represent the potential solution is afictitious density function σ called a single layer potential,Φ(~x) =∫Γσ(s)G(~x,~y(s))dS,~x ∈Ω. (6.5)The kernel G is the fundamental solution to Laplace’s equation in three dimensions,G(~x,~y) = 14pi|~x−~y| .In the limit as ~x approaches the boundary, we have the following combinedfirst/second kind integral equation for σ given the respectively Dirichlet/Neumannboundary conditions on Γ:Φ(~x f ) =∫Γσ(~y(k))G(~x f ,~y(k))dk, (6.6)∂∂nΦ(~xs) =σ(~xs)2+∂∂ns∫Γσ(~y(k))G(~xs,~y(k))dk. (6.7)Once σ is calculated, we can evaluate the gradient of Φ on the surface bytaking the gradient of (6.5) and replacing the normal component with (6.7). Thisprocess of computing the gradient is straightforward and lends itself to our numer-ics. Attempting a similar procedure with the direct formulation leads to hyper-singular integrals, which are far more difficult to handle numerically especially oncoarse triangle meshes. Additionally, the indirect formulation sidesteps a numer-ical oddity of the direct formulation, where the discretization treats the Neumannand Dirichlet data with the same approximation basis.A primary characteristic of our formulation and those in the engineering lit-erature is that they involve solving a Fredholm integral equation of the first kind(6.6). Mathematically, this implies that finding σ may be somewhat ill-conditionedsince our kernel is weakly singular (roughly, weakly singular means the kernel issingular, but integrable in L2 ) and therefore the integral operator is compact. Inpractice, both we and the engineering literature find meaningful results from solv-ing such a system, though care must be taken when handling the singular integralsto avoid unsolvable systems, and even then some degree of high-frequency artifactsshould be expected and regularized or smoothed away. Regardless, handling both49free surface and solid boundaries with straightforward integral equations likely ne-cessitates solving a first kind equation if we want to avoid solving systems withhyper-singular kernels (again roughly, hyper-singular means the kernel singularityis not L2 integrable and that the integral operator must be considered in a specialmanner).Advantages of the boundary integral formulation include:• Resolving everything on the surface dimensionally reduces our problem.• Integral equations have some surprising numerical characteristics that furtherincrease this efficiency• The infinite depth of the ocean is implicitly satisfied by the mathematicalformulation.We will further elaborate on the numerical advantages in the following section.6.3 Numerical Method for Boundary Integral WavesTo numerically solve the wave equations we take the following pseudo-algorithm:1. Solve the BIEs (6.6) and (6.7) for σ .2. Calculate ∇Φ on Γ.3. Advance Γ in time using equation (6.4) .4. Advance Φ(~x f ) in time using equation (6.2).6.3.1 Step 1: Solving our BIEs by CollocationThe first stage in a time step is constructing the potential solution satisfying ourboundary conditions. We approximate Γ using a triangular mesh comprised of tri-angles Tj and vertices ~xi. We will refer to a per-triangle enumeration of verticesas ~v j1,~v j2,~v j3,∈ Tj. A collocation scheme simply requires that the integral equa-tions are satisfied at discrete locations: we choose the mesh vertices. This gives thefollowing semi-discrete system:50Φ(~xi f ) =∫Γσ(~y(k))G(~xi f ,~y(k))dk, (6.8)∂∂nΦ(~xis) =σ(~xis)2+∂∂n∫Γσ(~y(k))G(~xis,~y(k))dk. (6.9)To fully discretize the system requires numerical quadratures for the integrals. Forequation (6.9), the theoretical properties of a second kind integral equation meanswe can expect numerical stability and well-conditioning with standard quadraturethat smooths the singularity.In contrast, the first kind integral equation (6.8) poses immediate challenges toeffective numerical quadrature. The kernel is integrable with smoothness assump-tions on the geometry and σ , but unfortunately any discrete approximation willhave to deal with the singularity of the kernel. In addition, theoretically the inte-gral operator is compact and therefore inverting it is ill-conditioned. However, solong as appropriate resolution, regularization or filtering is used, we can still attainmeaningful results, following in the footsteps of prior numerical work [25, 33, 78].We will attempt to provide intuition for the issue.It is readily apparent that a kernel whose action smooths the argument σ willresult in a first kind equation for which the inverse is ill-posed, since very dis-tinct rough inputs could be smoothed to virtually identical smooth outputs. Thus,one expects that treating the singularity by a smooth approximation will resultin an extremely ill-conditioned, if not unsolvable, linear equation, producing ar-bitrary high-frequency noise even if the low-frequency smooth components arewell-determined. Numerically solving such a system, one would expect to see in-creasing high-frequency noise as the resolution increases as explained by Groetsch[34]. However, in regards to solving the weakly singular integral equation of thefirst kind such as ours, Atkinson [5] states that they can be mathematically well-behaved and gives examples of numerical treatment that result in a stable numericalsystem with accurate solution to the two-dimensional analog of our equation.It is not apparent how Xue et al. handle integrating the singularity in theirmethod, though they mention seeing fine scale noise that they accordingly filter;Fochesato et al. simply state that they use “recursive subdivision”. We accurately51integrate the singularity similar to Grilli et al. by first applying an integral transfor-mation to remove the singularity, then using Gaussian quadrature. This process iscommonly known as Duffy quadrature [20]. In our case, we found it resulted in awell-behaved linear system that gave results consistent with our low-order geom-etry approximation. In contrast, neglecting the singularity and using standard loworder quadrature for the first kind equation resulted in our linear solver failing toconverge.Figure 6.3: Ocean Waves initialized with Phillips spectrum [71].Numerical QuadratureOur first step in implementing our quadrature is to approximate σ by discrete val-ues σi at each mesh vertex, linearly interpolated across triangular faces.We calculate the equations (6.8) and (6.9) at a mesh vertex ~xi by integratingover all triangles. For those triangles which don’t include ~xi we use a simple 3-point quadrature rule,∫Tjσ(~y(s))G(~xi,~y(s))ds≈ 13[∑k=1,2,3σ(v jk)4pi|~xi−~v jk|)]A(Tj) xi /∈ Tj, (6.10)and52∫Tjσ(~y(k))∂∂nG(~xi,~y(k))dk≈ 13[∑k=1,2,3−σ(v jk)~nis · (~xis−~v jk)4pi|~xis−~v jk|3 )]A(Tj) xi /∈ Tj, (6.11)where A(Tj) is the area of triangle j. When vertex ~xi ∈ Tj, and we are calculatingequation (6.9) we adapt the 3-point quadrature from equation (6.11) by setting σito zero at ~xi =~vi j, giving an O(h2) approximation for characteristic edge lengthh. This low order approximation is acceptable since the diagonal elements of theresulting matrix are dominated in the limit by the jump condition in (6.9). Whenvertex ~xi ∈ Tj and equation (6.8) is being calculated, we instead use Duffy quadra-ture, described below, to handle the singularity.Duffy QuadratureDuffy quadrature uses a simple variable transformation to completely remove thesingularity from the integrand. For the sake of simplicity, we describe Duffyquadrature on a reference triangle Tr in the x-y plane, vr1 = (0,0,0), vr2 = (1,1,0),and vr3 = (0,0,1), where the singularity is at vr1. All other triangles can be affinelytransformed to this reference triangle. On the reference triangle we can write theright hand side of (6.6) as the 2D integral14pi∫ 10∫ x0σ(x,y)√x2+ y2dxdy. (6.12)With the transformation, u = y/x, we arrive at14pi∫ 10∫ 10σ(x,xu)√x2+(xu)2xdxdu=14pi∫ 10∫ 10σ(x,xu)√1+(u)2dxdu. (6.13)53The singularity in equation (6.12) no longer exists in (6.13). We approximate (6.13)with a nine-point tensor Gaussian quadrature on the square, using the linearly-interpolated σ . The linear transformation into the reference triangle of courseneeds to be taken into account in general.Linear SolverThe result of applying quadrature to equations (6.8) and (6.9) is a dense linearsystem. Actually constructing the dense matrix for a direct solve would be unac-ceptably expensive. We instead look to iterative solvers that don’t require the wholematrix but merely compute several matrix-vector products, choosing BiCGSTAB[72] in our implementation. Moreover, we can further accelerate matrix-vectormultiplication with parallelism and the Fast Multipole Method, described in sec-tion 6.4. For smooth free-surface-only geometry we found BiCGSTAB solved thesystems to 10−5 of the initial residual norm in less than 30 iterations. With solidgeometry included, it typically took 50–60 iterations. However, when the meshbecomes overly distorted with sliver triangles, the iteration count can rise past 100,and if the simulation becomes unstable and the geometry pathological, the solvercan fail to converge. We thus employ the El Topo library to maintain a well-formedmesh by restricting maximum and minimum edge lengths, and further stabilize thenumerical method by adding some viscosity described below.6.3.2 Step 2: Computing ∇ΦThe normal derivative at the surface is given by equation (6.7) and can be calculatedby applying the same process describe for computing (6.9). The derivatives in thetangential directions~t we compute by applying the gradient directly to equation(6.5). There is no jump in the tangential derivatives as they approach the boundary,yielding~ti ·∇Φ(~xiΓ) =~t ·∇∫Γσ(~y(k))G(~xiΓ,~y(k))dk. (6.14)We compute this in the same low-order fashion as before,54~ti ·∇∫Tjσ(~y(k))G(~xi,~y(k))dk≈ 13[∑k=1,2,3−σ(v jk)~ti · (~xiΓ−~v jk)4pi|~xiΓ−~v jk|3 )]A(Tj) xi /∈ Tj. (6.15)Again we neglect the singular components when ~xi ∈ Tj by setting σi to zero there.We compute the tangential derivative for two perpendicular tangential componentsthen combine them with the already computed normal component to construct thegradient of the potential on the boundary. The engineering literature computesthe tangential gradients by making a local high-order approximation to the surfacefrom which they can approximate the tangential derivatives using the necessarygeometric information and values of Φ. See especially the method by Grilli et al.[33]. Our method is a far simpler way of constructing the tangential derivative andmore robust for unstructured triangle meshes where local high-order approximationto the surface would be a challenge.6.3.3 Steps 3 and 4: Time IntegrationWe update the surface and potential at time t = ∆t · (n+1) using symplectic Euler:~xn+1i =~xni +∆t∇Φni (6.16)Φn+1i =Φni +∆t[p(~xn+1i )−gz(~xn+1i )]. (6.17)We further stabilize the system by including (in essence) some artificial viscosity,adaptively smoothing the mesh in the tangential plane. This is done by computinga weighted average of each vertex’s position with its m neighbors,~˜xi = α~xi+βm∑k=1~xk, (6.18)β = (1−α)/m, (6.19)55and then moving the vertex to the projection of the weighted average on thetangent plane. We found that smoothing in the tangential plane resulted in wellformed triangles without noticeable energy loss due to the altitude erosion that onewould see with regular smoothing. On the solid surface the boundary values of Φis not propagated in the model and the fluid geometry can be changed per time-step without affecting the solution. Here we tangentially smooth quite stronglywith values of α = 0.05. On the free surface, we generally smooth with α valuesaround .97 for a time-step of dt = .01, computing four simulation time-steps perframe. At the edge of our simulation, we smooth Φ and the mesh in all planes toavoid problems when large waves hit the mesh edge. At the edge points we smoothonly in the y-coordinate, but still allow motion due to the potential gradient.6.3.4 Solid Boundaries and Mesh AdaptivityTo determine whether a fluid point is a solid or free boundary, we check auxiliarymeshes denoting solid objects for collision with the fluid mesh. Potential collisionsare handled by the El Topo mesh library. A fluid point that is in an adjustableproximity to a solid becomes considered part of the solid boundary and is movedto a pre-defined distance from the solid. Solid normal velocities are then applied tothis vertex during time-stepping and it is considered to have the requisite Neumannboundary condition in the numerical solver.For deep ocean waves, we found that the natural tangential motion of the wavedynamics caused the compression of the surface mesh near areas of high curvaturesuch as the peaks of the waves. This gives an added efficiency to the method,imparting an effective resolution many times higher then the original mesh size.While this adaptivity can sometimes be reversed due to unnatural pressure forcing,it is consistent for areas where the natural wave motion is allowed to occur.6.4 FMM and GPU accelerationBoundary integral equations reduce the geometry and thus memory requirementsof wave simulation, from a 3D grid to just the 2D surface mesh, and exhibit niceconvergence properties with using iterative solvers. However, the O(N2) com-plexity of multiplying with the (conceptually) dense matrix means we’ve essen-56Figure 6.4: Sequence of images showing waves with longer wavelengthovertaking those with smaller, a characteristic of deep water dispersion.tially moved complexity from space into time. Fortunately, we can reduce the highcompute cost of solving boundary integral methods. We employ two different ap-57Figure 6.5: Top down view of the computational mesh showing natural adap-tivity induced by Lagrangian motion.Dims Verts FMM time GPU time80x80 6400 6.6 6.0E-2160x160 25600 30 6.3E-1320x320 102400 96 9.0640x640 409600 3.5E2 1.3E2Table 6.1: Seconds per BiCGSTAB iteration using FMM or GPU accelera-tion for mesh dimensions (Dims) and total number of vertices (Verts)proaches, using the FMM and also GPU acceleration to more efficiently computethe wave simulation. In both cases, instead of computing integrals by summing pertriangle as described, we first compute all the quadrature weights accumulated by avertex from its neighboring triangles, and then compute the integrals by summingover each vertex. Computing the quadratures in this fashion becomes an N-bodyproblem involving the fundamental solution for Laplace’s equation and its gradient.This can be done in O(N) time using the FMM. In addition, the N-body problem isembarrassingly parallel and is straightforward to implement in OpenCL to exploitthe massive parallelism of modern GPUs.58Figure 6.6: Logarithmic plot of table 6.1 showing the number of verticesversus time per iteration asymptotics of the FMM and the GPU.The FMM we used ran on a single CPU thread on a 3.2Ghz Intel i5-3470 at 3-4digits of accuracy; though it showed linear-time asymptotics, as seen in table 6.1and Figure 6.6, it was still slower on our examples then the direct GPU-acceleratedversion using an AMD Radeon 7970. We are currently working on a GPU versionof the FMM to combine these two acceleration techniques for added computationalefficiency.An additional advantage of our BIE formulation is that it allows us to easilyuse our FMM, implemented according to that described in [15]. The formula-tions in the engineering literature require dipole charges rather then single pointsources, which necessitates non-trivial modifications to the numerical machineryof the usual FMM.59Figure 6.7: Waves incident upon a near-vertical wall.Figure 6.8: Moving balls with wake.6.5 ResultsOur results demonstrate deep ocean wave phenomena. Figure 6.1 shows the steeppeaks and shallow troughs expected for ocean waves. Figure 6.3 shows a largescale simulation of waves initialized with a Phillips spectrum [71]. The dispersionrelationship for deep water waves imparts higher wave speed to waves of longerwavelength. Figure 6.4 shows this happening in our simulation. While this visuallyapparent dispersion relationship is inherent in FFT simulations and evident in ours,60it is not captured by shallow water simulations. Figure 6.9 shows a comparisonwith a height-field waves simulation using FFT’s. Figure 6.7 shows waves crash-ing on a near vertical wall, demonstrating solid boundary interactions. Figure 6.8shows two balls moving on the water surface, demonstrating moving solid bound-ary interactions. Figure 6.5 demonstrates the natural adaptivity of our method,allowing a regular grid to distort naturally, compressing the grid at the sharp peaksof the waves and decompressing it in the smoother troughs.Figure 6.9: Comparison of our waves (top) with a Tessendorf/FFT height-field wave simulation (bottom) using the same initial conditions. Color-ing represents the potential values on the surface .616.6 Future WorkOur BIE method is a novel method for deep ocean wave animation, and thus thereare ample avenues for future work:• The greatest immediate barrier to effectiveness of the method is its relativelyhigh compute cost. As we’ve shown, the FMM and GPU computing canameliorate that. Combining the two could bring us closer to large-scale,low-memory ocean simulations with dynamic solid boundaries.• The method is dependent on good mesh-tracking software. While El Topowas well suited to the current method, phenomena such as overturning waveswill need meshing algorithms that handle topology changes robustly and im-mediately, without leaving thin interfaces of air between fluids as El Topo iswont to do.• Stable ways to incorporate the non-linear term would expand the use casesof the method. While the non-linear term is not necessary for general oceansimulation, specific phenomena such as overturning waves and wave-breakingwould require implementing this term and provide interesting problems intheir own rights.• Tight integration with Tessendorf waves to continue ocean simulations outto the horizon would be a straightforward next step. It is simple to derive thesurface potential values from a FFT based wave simulation; one could easilyuse this information to drive the boundaries of our method. This would allowseamless interfaces between the more dynamic and interactive BIE methodand a periodically extended FFT-based procedural ocean.• Non-reflecting conditions on the edge of the domain where the free surfacecontinues.• Two-way coupling with buoyant objects would allow more dynamic oceanscenes. This would require developing a buoyancy model compatible withpotential flow.62Chapter 7Compact Iso-SurfaceRepresentation and TemporalCompression for FluidPhenomena7.1 SynopsisWe propose a novel method of compressing a fluid effect for real-time playback byusing a compact mathematical representation of the spatio-temporal fluid surface.To create the surface representation we use as input a set of fluid meshes fromstandard techniques along with the simulation’s surface velocity to construct a spa-tially adaptive and temporally coherent Lagrangian least-squares representation ofthe surface. We then compress the Lagrangian point data using a technique calledFourier extensions for further compression gains. The resulting surface is easilydecompressed and amenable to being evaluated in parallel. We demonstrate real-time decompression and meshing of the surface using a dual-contouring methodthat efficiently uses the decompressed particle data and least-squares representa-tion to create a view dependent triangulation.637.2 IntroductionThe output meshes of modern fluid simulations are highly detailed but temporallyincoherent due to standard surfacing techniques for volumetric methods or meshmerging and refinement for ones employing explicit mesh tracking. The combi-nation of large mesh size with temporally incoherent triangulation decreases theutility of these meshes for real-time use cases such as games or virtual reality dueto storage and bandwidth issues; naively, connectivity and point data must be re-tained for every frame. While storage and bandwidth are of concern to real-timeapplications, adaptive reconstruction or playback of data with varying levels ofdetail is also a primary concern. Varying the level of detail of triangulated meshdata requires creating new lower resolution meshes which create even more data,or unwieldy runtime decimation schemes. For standard non-temporally-coherenttriangulated data, the standard approach is to apply a spatial mesh compression foreach sequential mesh per time-step. Some methods include connectivity compres-sion; if the connectivity changes each time step, it must be re-compressed eachtime.In terms of storage and bandwidth of fluid animation, these triangle basedmethods are inherently inefficient since the fluid itself is a product of continuousphysical forces and exhibits significant temporal coherency. The temporal incoher-ence of the triangular output indicates that the triangle connectivity itself is of littleimportance to the utility of the surface. The object of primary importance is actu-ally the geometry of the implicit surface — the triangular mesh is just a discreterepresentation necessary for rendering.We propose a novel method of compressing temporally coherent time depen-dent iso-surfaces for real-time playback. We create a compact mathematical repre-sentation of the spatio-temporal surface using a Lagrangian set of surface incidentpoints and their normals. The representation consists of a least squares approx-imation consisting of blended second order polynomials which are calculated byexplicitly approximating the surface curvature at each point, dependent on the ge-ometry of that point’s neighbouring surface points. We construct this representa-tion using triangulated mesh output from standard simulation techniques, with theaddition of the values of the simulation’s per-vertex surface velocity. The repre-64sentation is spatially adaptive and temporally coherent. The spatial least squaresapproximation gives a super-quadratic lossy spatial compression while the coher-ence allows us to apply additional geometric temporal compression by representingthe time data with a high-order interpolation called Fourier extensions.The demands of real-time effects requires at least an ability to reduce compu-tational and rendering costs as the effect recedes in proximity to the viewer. Asstated, this is difficult to do efficiently using direct triangular mesh formats. Ourmethod enables creating the triangular geometry in real-time, allowing the finalgeometry to be directed by the view point. We demonstrate meshing the surface inreal-time using a dual-contouring method that efficiently uses the particle data andleast-squares representation to create a view dependent triangulation.Significant work has been done on triangle mesh compression. We refer thereaders to the review paper by Maglo et al. [53]. Most of these schemes don’tcompare directly to our work in that they are primarily concerned with triangulatedgeometry, rather than the implicit surface itself. For temporal mesh compression,in general these techniques require temporal coherency in the mesh connectivity;as stated above, this is unnecessary for fluid effects, and certainly unreasonable toexpect for compressing output of standard fluid techniques.Of substantial note is the method of Karni et al. [40] where they define a spec-tral implicit surface approximation by segmenting the triangular mesh and usingFourier patches as surface representations; this technique is still used to computetriangle geometry and connectivity compression. One could utilize the spectralapproximation of Karni et al. in the same manner as we do by using the spectralsurface approximation as the compressed stored representation of the surface. Ourspatial representation, though of lower order approximating power, maintains a La-grangian representation of the surface that is amenable to temporal compression.This would be difficult to do with Karni et al.’s segmenting Fourier-based tech-niques; the segmenting would be difficult to maintain as a temporal componentunder topology changes. Additionally, our experience using spectrally accuratemethods (radial basis functions) for the spatial representation leads to large densesystems with larger temporal variance in high frequency coefficients, making themless suitable for temporal compression. We believe that in general higher orderspatial approximation techniques are unsuitable to temporal compression for this65reason.We would also make reference to the real-time smoke lighting algorithm cre-ated by Zhou et al. [82]. By using significant precomputation, they were able torepresent initial volume data in a compressed fashion that also enabled run-timebenefits; real-time volumetric lighting in their case. In the same manner, we viewour “compression scheme” not only as a method to reduce data size, but also asa compact representation which enables fast and adaptive parallel triangulation ofthe underlying surface at run-time.7.3 Lagrangian Surface RepresentationWe base our implicit surface representation on a surface distribution of point posi-tions~xi and their normals~ni that move with the material motion of the surface. Wewill refer to these points as least-squares points (LSQ points). Using this data wefirst construct a local approximation to the surface at each LSQ point using neigh-bouring points and then blend that per point approximations to fully construct thesurface. For practical use, we consider the triangulated data from a standard fluidsimulation with normal and velocity data at the vertices as the representation ofthe implicit surface that we are attempting to approximate. We will refer to thissequence of triangle mesh input as submeshes.7.3.1 Pointwise ReconstructionFor each LSQ point position ~xi we construct an implicit function Γi(~x) for whichthe zero contour approximates the surface in the points vicinity. This pointwisereconstruction is defined asΓi(~x) = d+~ni ·~ri+aτ2i1+bτi1τi2+ cτ2i2 (7.1)where a,b,c,d are computed geometric parameters, ~ri is ~x−~xi and τ1 and τ2 arethe projection of ~r on unit tangent vectors ~τi1,2 that are orthogonal to each other66and the normal~ni at~xi.τi1 =~ri ·~τi1τi2 =~ri ·~τi2The first two terms of our reconstruction (7.1) constitute a planar approximationto the surface; in order for the zero contour to be at the point’s position, d is zeroby construction. Along the unit vector from the LSQ point position we have anapproximation to the signed distance function. The quadratic terms of the functionallow us to account for curvature. If we set Γ to zero and solve for the normalcomponent of~ri , then (7.1) constitutes a 2nd order Taylor series approximation tothe surface, assuming the coefficients a,b,c are the appropriate second derivatives.Since we don’t have a functional representation of the surface to calculate deriva-tives with, we resort instead to minimizing a least squares problem for a,b,c thatattempts to minimize the value of Γi(~x) at neighbouring LSQ points. Due to ourapproximation of the 2nd order coefficients, we don’t actually expect to achieve3rd order accuracy, however, we do empirically find that this reconstruction is atleast super-quadratic.Given a neighbouring set of LSQ points ~x j, we construct a weighted (usuallyoverdetermined) linear system with a, b, c, and coverage radius δcΓi(~x j) = 0 =(~ni ·~ri j +aτ2i j1+bτi j1τi j2+ cτ2i j2)W (2|~ri|/δc) (7.2)W (s) =(1− s2)2 0 < s < 10 1 < sWe solve the normal equations for the overdetermined system (7.2). Constructing(7.2) requires at least three neighbours. Reliably approximating the surface curva-ture by constructing values for a,b,c depends on determining a sufficient neighbourdistribution. We will cover this in the next section.7.3.2 Global Distribution and ReconstructionWe create a distribution of LSQ points across a surface with an average charac-teristic geodesic distance between points called the coverage radius (described in67section 7.3.4). Reconstructing the global surface becomes a question of interpo-lating between the pointwise reconstructions. An interpolation method needs torespect that as one approaches a LSQ point on the surface, the reconstruction ofthe local pointwise function must be the limiting one, since we know it is exactlyplaced on the surface. A standard way to do this would be to use an inverse dis-tance weighting, similar to Shepard’s [68]. We found that this method was proneto rounding errors close to the points.We instead use a method reminiscent of moving least squares surface con-struction as described by Levin [49], though our implementation is somewhat sim-pler. We interpolate the surface by combining the pointwise reconstructions as aweighted average. We use those pointwise reconstructions contained in a variablesearch radius set usually to two times the coverage radius around the evaluationpoint. The global implicit function Ω is computed asΩ(~x) =1∑iW (|~ri|/δ )∑iW (|~ri|/δ )Γi(~x) (7.3)W (s) =(1− s2)2 0 < s < 10 1 < s (7.4)δ is a scaling factor dependent on the two closest neighbors to the evaluation point~r1,~r2 with |~r1|< |~r2| .δ = |~r1|+0.9|~r2|. (7.5)We bound δ at two times the coverage radius. This dynamic scaling factor forcesthe interpolation to converge to the nearest pointwise reconstruction as the evalua-tion point moves towards it. The influence of the rest of the points~xr are diminishedas the scale factor drives the value of |~rr/δ | towards and above one, which gives asmaller and eventually zero W . The result is a global interpolating function with-out the instabilities of singular weighting kernels that converges to the pointwiseapproximations as we approach them on the surface.Our global interpolation scheme allows us to easily define the global surfacewhen point densities correlate with the coverage radius. The dynamic scaling fac-68tor enables smooth approximation even when Lagrangian movement causes theparticle position distribution to be less uniform, especially where point proximitiesbecome small.7.3.3 Sharp FeaturesBy using a polynomial approximation of the implicit surface, we can accuratelyand efficiently approximate smooth surfaces. However, sharp features are not ef-ficiently captured by such an approximation. Fleishman et al. [24] used statisticalmethods to derive normal data from point cloud distributions which informed asharp feature least squares implicit construction from scattered point data. In ourreconstruction, we have a-priori normal information from the original mesh, andcan thus detect and resolve sharp features in the input data. We augment our La-grangian least squares reconstruction to utilize this normal data to preserve sharpedges in the reconstruction.Our pointwise surface construction (7.1) is valid as a local approximation to asmooth surface. On or near a sharp feature, this approximation becomes less valid;for example, on convex edges a polynomial approximation leads to a local flatten-ing and lateral bulging along the feature due to the radial nature of the polynomialapproximation as seen in figure 7.1. To avoid this, we do not place LSQ points onmesh vertices with adjacent face normals whose relative angles exceed a control-lable value. To maintain a smooth approximation during coefficient calculation foran LSQ point when close to a sharp feature, we don’t include neighbouring LSQpoints whose normals are far from parallel when solving the coefficients.Then, when reconstructing the LSQ surface at an evaluation point, we detectwhen neighbouring LSQ point normals are not close to parallel by calculating howfar the magnitude of the average of neighbouring normals is from one. Insteadof summing the neighbouring least squares reconstructions as usual, we insteadconstruct the implicit value by summing positive and negative values from neigh-bouring LSQ points separately. For convex features, we decrease the negative con-tributions by a multiplying by a user-defined power of the normal average magni-tude. This provides a local sharper ”boolean” reconstruction to the surface, wherethe two planar approximations converging to a sharp feature are treated similar to69space partitioning planes, as seen in figure 7.1.Figure 7.1: Resolving sharp features on a cube. On the left we are using ourweighted normal averaging. On the right, no sharp feature resolving isapplied.7.3.4 Particle DistributionWe have described how we reconstruct an implicit surface given a valid distribu-tion of points. In this section we will detail how we construct this distribution fromthe initial triangular fluid meshes, which we call submeshes. We start by makingan initial distribution by randomly choosing points off the submesh at random, anddiscarding those that are within the coverage radius distance to another point. Thisdistribution doesn’t take into account the geometry of the fluid surface, which maybe thin or close to overlapping; We therefore modify the distribution so that it isdefined by a quasi-geodesic distance (QGD). The QGD is defined by finding theshortest combination of connected edges between two points, and then summingtheir lengths. This QGD Voronoi diagram is calculated on the submesh by march-ing along the edges from the already chosen points. Each step of the marchingalgorithm across an edge requires calculating the distance moved along that edgefrom the previous point. The new point’s QGD is then calculated as the previouspoint’s plus the edge distance. The edge steps are queued up and processed inorder of least total QGD. When a step encounters a submesh vertex that alreadyhas a QGD value, it terminates when that distance is less then the newly calcu-lated one. Once this initial weighted graph distance is calculated, we then start70adding submesh points to the LSQ point distribution, starting with those with thehighest QGD from an initial point. Points are added greedily at submesh verticeswith largest QGD. Each time we add a point we update the QGD map accordingly,by marching outwards from the new point. The process is halted when the high-est QGD distance falls below the coverage radius. Figure 7.5 shows successivelydenser distributions on a smoothed cube initial geometry, shown in Figure 7.4.7.3.5 Spatial CompressionUsing a super-quadratic global approximating function allows us to apply spatialcompression when compared to the O(h) approximating powers of triangulatingthe surface. Additionally, our reconstruction scheme needs only point and normaldata. The higher order coefficients are calculated explicitly and locally at run-time,constituting additional compression. The additional advantages of calculating thehigher order data at runtime is that the representation we store is the lower orderposition and normal values; these values are expected to be temporally coherent,and allow us to apply compression in time.7.3.6 Lagrangian TransportOnce we have constructed the LSQ point distribution corresponding to the surfacerepresentation of the initial frame’s submesh, we then determine velocity valuesof those points via barycentric interpolation from the submesh vertices and thenadvect them with a simple Euler step using the frame’s timestep. The new advectedpoints are then “snapped” to the next frame’s submesh by calculating the closesttriangle and projecting the point onto it. We then repeat the process described inthe previous section, new points are added by calculating a QGD function thatdetects areas lacking point coverage. To combat oversampling when the surfacevelocities naturally bring particles closer together, we assign a finite lifetime foreach particle. Any resulting gap in the particle coverage is allowed to be filled bythe particle distribution mechanics we have already described.717.4 Temporal CompressionThe points of our Lagrangian surface representation follow the material flow gen-erated by the fluid simulation. In the absence of boundaries, this flow is a resultof continuous processes and is therefore expected to exhibit temporal smoothness.Our goal is to exploit this smoothness by compressing the time dependent data byrepresenting it by a spectral interpolation basis. This is not as straightforward as itsounds. Each particle’s data has a finite lifetime defined by equispaced time-steps.Both simple Fourier and polynomial interpolation of this kind of data experienceartifacts because of periodic discontinuities or the Runge phenomenon, respec-tively. We use instead a recent advancement in mathematics called Fourier exten-sions which provide us with a spectral interpolation of our equispaced time datawithout the problems of standard spectral interpolation methods. Fourier exten-sions are simply the least-squares projection of interpolated data on a sub-intervalof a Fourier interpolant. In our case, we scale our time dependent data f (ti) fori= 0, . . . , j onto an interval of [0..12 ] and then calculate its projection onto a trigono-metric Fourier series consisting of a0+∑nk=1 Akeik2pit with a0 and Ak being real andcomplex coefficients, respectively. We use a simple collocation scheme to con-struct the following system of equations:a0+n∑k=1Akeik2piti = f (ti) ti =i2 jfor i = 0 . . . j (7.6)We achieve compression by looking for n < j, and since our Fourier series isn’tnecessarily orthogonal on this interval, we resort to solving this equation by usinga standard least-squares algorithm from a readily available numerical software li-brary. We initially start with small n and increase it until a user-defined bound onthe l1 error is satisfied by the reconstructed function. The Fourier series now accu-rately approximates our time data on [0..12 ] and is still able to satisfy periodicity onthe completed [0..1] interval.Astute readers will realize that we are calculating the L2 projection of our time-data on a sub-interval of the usual domain over which the Fourier series is orthog-onal and that on this sub-interval the underlying basis functions are no longer or-thogonal. This means that our coefficient matrix in equation (7.6) will tend to be72ill-conditioned, this necessitates solving it using standard SVD-based least-squaresmethods. In order to be efficient, we allow our off the shelf solver to ignore in-significant singular values in its computations of the projection. Adcock et al. [2]showed that this method of computing extensions is actually stable and that theextension converges at least super-algebraically.The final per-particle compressed data is the set of Fourier coefficients approx-imating the time-dependent data for each vector component of the position andnormal. This floating point data is stored in sequence by the start frame of eachparticle. Some additional meta-data is necessary to reconstruct the particles. Westore three integer values per point, the beginning frame, the lifetime in frames ofthe particle and the index for its relevant Fourier data in the floating point array.Each of the six Fourier series can be adaptively sized, and we keep an additionalinteger for each to indicate the length of the Fourier series.7.5 Decompression and Re-MeshingThe LSQ points and their normal data are decompressed simply by evaluating theappropriate Fourier data at the correct time. The local curvature values are thencalculated. At this point, the global reconstruction can be used to evaluate the sur-face. Most applications, however, require a triangulated surface, rather then a levelset function. For this purpose, any triangulation technique, such as the widely usedmarching cubes algorithm [51] could be used with our surface representation togenerate the required data. Our Lagrangian point data set itself is by design colo-cated on the implicit surface. Similar to Wyvill et al. [77] we utilize the underlyingstructure of our surface representation to efficiently mesh it. We describe a dualcontouring technique that efficiently uses the point positions which are incident onthe surface. The initial surface intersections are calculated at point positions, andthen subsequent triangles are calculated by marching outwards from these initialcells, again similar to Wyvill et al. [77]. The method is easily amenable to paral-lelization, a critical requirement of any real-time triangle construction, consideringthe parallel architecture of today’s graphics hardware.We begin by defining a grid system with cell length determined by the desiredtriangulation resolution. We then instantiate each grid cell that contains a decom-73pressed LSQ point. The implicit surface is evaluated on the grid corners of thesegrid cells. Sign switching among any of these corners signifies a contour intersec-tion on the cell itself. By the nature of our reconstruction, each LSQ point mustlie on the surface. It can happen, however, that the LSQ point is interior to a cell,but that the surface only intersects along a single face. We detect this case andthen evaluate the neighbouring cube instead. If that cell also does not have a signchange, it is discarded, otherwise it is associated with the LSQ point.As in other dual contouring techniques, each grid cell with a sign change isassociated with an interior point. In the case of the initial grid cells, this point willbe the LSQ points. In the case that we have multiple LSQ points associated with asingle cell, the point nearest to the cell center is chosen.From the evaluated values on the cell corners, we determine the neighbour-ing cells intersected by the implicit function and also calculate edge intersectionpoints using a linear approximation along the edges. We then repeat the process ofinstantiating and evaluating neighbouring cells, with the added step that for eachnew cell we calculate a new interior surface point based on linearly interpolat-ing the edge intersection points. The marching process continues as each newlycreated cell also queues non-processed intersected neighbour cells, until either nonon-processed intersected neighbours are available (if the contour is manifold) ora user defined limit on marching iterations is reached. The triangulation is finallycreated by creating four element triangle fans around each edge, with the edge in-tersection point being the pivot vertex, and the four cell points adjacent to the edgebeing the vertices on the outer cycle.This method of triangulation is by design low order; we could easily implementmore accurate edge and cell intersection points, for example by following Ju et al.[39]. However, it does exemplify some of the advantages of our implicit surfacerepresentation. Information about the surface is given immediately by the positionof the decompressed points. Our marching technique is easily modified to run inparallel. For real-time application, if we set the grid cell length to be roughly thatof half the coverage radius, we are then guaranteed marching will succeed in just acouple iterations. At this point, using the already given cell points offers additionalcomputation saving.Creating the triangles at run-time offers significant control over the level-of-74detail being rendered. The first and simplest way is to simply control the size of thetriangulation based on camera proximity. In addition, we could control the pointdistribution, simply by dynamically choosing a subset of the points with a largercoverage radius. We could also implement our triangulation to be adaptive, as in Juet al. [39] so that for large complicated surfaces, we could apply our level-of-detailtechniques gracefully on separate portions of the same mesh. Our method wouldalso be amenable to screen space or frustum conforming techniques, for example,that described by Mu¨ller et al. [56].7.6 ResultsFigure 7.2: The meshed torus, white points show the least squares positionsand extruding red lines showing the analytic normal direction.Our compression technique is a discrete representation of a time-dependentimplicit surface. A triangulated mesh sequence, along with mesh compression75techniques, can be considered similarly, leading to the possibility of comparisonof our method to triangular mesh compression techniques. However, comparisonwith mesh compression techniques isn’t immediately straightforward since mostresults for these methods consider the vertex distance from the original mesh tri-angulation, rather than the implicit surface which it approximates. This meansthat traditional mesh compression includes an inherent first order approximationerror, whereas ours is super-quadratic. We demonstrate this approximation powerby figures 7.3 and 7.2, which shows the super-quadratic convergence of our leastsquares approximation to the analytic representation of the torus S. To constructthis approximation, we first calculate a dense triangulation of the torus using adual contouring technique that employed a numerical root solver to place each ofits triangulation vertices on the torus’ surface within rounding error. The normalswere then produced by computing the analytic gradient. From this approximation,we selected a subset of the mesh vertices and normals as our least squares-points,using a rough Poisson distribution which maintained a point distance h. To com-pute the error of the least square implementation, we calculated the distance of thereconstructed surface S¯ from each of the vertices of the dense initial triangulation,again by root finding. By taking the greatest distance, we arrive at an approximateL∞ norm for the triangulation. Given that we know the curvature estimates of thetorus and the least-squares approximation, we can be certain of our error estimate,rather then resorting to more thorough Hadamard distance measures. As figure 7.3shows, our least squares approximation achieves super-quadratic convergence inapproximating the analytic torus.Table 7.1 lists the compression rates for 240 frames of a smoothed cube shownin Figure 7.4 undergoing rotation and translation. The smooth cube consisted of1538 vertices and 3071 faces and exhibits a smooth surface with varying curvature.The LSQ point distribution corresponds to the coarsest example in Figure 7.5. Theamount of temporal compression is determined by setting the required temporal er-ror (RTE). In our implementation, the Fourier extensions will reduce the number ofinterpolating basis functions until this error is surpassed for compressing the timedependent values of positions and normals of the LSQ points. The average errorand maximum error show the normal distance from the original mesh vertices tothe approximating surface, the smoothed cube’s length dimensions were approx-76Figure 7.3: Convergence of the least squares representation on the surface ofan analytic torus. The log plot shows the error convergence with a linearfit for the approximate asymptotic order.imately 2 units. Total, Spatial and No Compression (TC,SC,NC) list the numberof words; i.e. the number of floats and integers, required to represent the mesh inits totally compressed form, only spatial compression using per Frame LSQ pointsand using the original per frame vertex and face data directly. The smooth cubeexample adaptively created triangles dependent on view proximity to the object.For 12,000 triangles from 100 LSQ points the replay ran at 50 frames per secondon an Intel Core I5.Table 7.2 gives the time required to compute the Lagrangian LSQ points dis-tribution (spatial compression) for all 250 frames of our test data, then computethe adaptive Fourier approximations of their time-dependent positions and normals(temporal compression). The decompression timings are the rough average of mil-liseconds used to compute the inverse Fourier transform in order to pass the LSQpoints into the real-time surfacing algorithm. Figure 7.1 shows the effects of re-solving sharp edges using the techniques we described in Section 7.3.3. In thiscase, we modified the negative component of the least squares summation by thefourth power of the averaged normal length.77RTE avg err max err TC SC NC.001 .0071 .063 54k 132k 3429k.01 .0073 .063 38k 132k 3429k.1 .014 .12 27k 132k 3429kTable 7.1: Spatial and temporal compression for a moving smoothed cube.CompressionSpatial .001 .01 .1298s 6.1s 4.1s 2.3sDecompression- ≈3.1ms ≈2.3ms ≈1.3msTable 7.2: Timings for spatial and temporal Compression and temporal De-compression.Corto .001 .01 .12218b 843b 614b 422bTable 7.3: Moving smoothed cube per-frame file sizes in bytes compared toCorto.Figure 7.3 compares the per-frame file size of our compressed data with thatproduced by the Corto mesh compression library [1] [63], compiled and run withdefault parameters. The per-frame values for the LSQSURF algorithm was calcu-lated as the total file size divided by the number of frames that it represented. Thefile size for each entry in Figure 7.3 corresponds to the binary data with RTE errordescribed in Figure 7.1. The smooth example we are approximating is advanta-geous to the LSQSURF algorithm. Exhaustive testing and comparison with differ-ent geometry types and deforming topologies would be interesting future work.78Figure 7.4: Smoothed cube test geometry.Figure 7.6 shows an interactive replay of a compressed mesh created by a com-mercial visual effects software. The single threaded decompression and meshingcreated 100,000 triangles from 1,500 LSQ nodes at 5 frames per second on an In-tel Core I5. For both of our replay timings, the majority of computation time wasspent in evaluating the surface during dual contouring.7.7 ConclusionWe have demonstrated that our Lagrangian least squares method of iso-surfacecompression and playback is easily parallelized, amenable to adaptivity and ableto capture sharp features. Further work could be done to increase the effectivenessof these aspects. Our method of distribution and tracking of LSQ points couldbe refined to place more points on geometry with sharper features. This could beinformed by error estimates of the intermediate distributions of points. Also, welimit the temporal compression by imposing a maximum number frame for the lifeof an LSQ points. It would be interesting to pursue a global optimisation strategythat had longer life-times in spatio-temporal areas of the animation sequence thatexhibited slow time variance. The sharp feature strategy we demonstrated could79Figure 7.5: LSQ points representations of a smoothed cube.Figure 7.6: LSQ points representation of a fluid simulation created fromcommercial software.80also be adapted easily to concave and thin features in the mesh. Finally, our dual-contouring method still requires a background grid, and the marching strategy isstill somewhat limited by its greedy implementation. For real-time playback, anyadvances in constructing the rendering mesh swiftly and efficiently would increaseits utility.81Chapter 8ConclusionsThe surface based animation techniques described in this thesis reduce both thetime and space complexity required to create detailed visual effects for smoke andwater. The combination of dimensional reduction due to boundary integral basedmodelling and numerics along with the modern algorithmic advancement of thefast-multipole method result in techniques that demonstrate efficient spatial andtemporal asymptotics when compared to equivalent volumetric based techniques.Adding efficient boundaries using a boundary integral representation describ-ing potential flow in chapter 5 provided a solution to a problem that had plaguedprevious vorticity based animation techniques.The boundary integral equation formulation for ocean wave dynamics in chap-ter 6 provided a fully dynamic surface-only animation technique that implicitlyincluded the infinite depth of the ocean, agreed with procedural techniques andprovided artist control with O(N) asymptotic complexity.Finally, in chapter 7, the Lagrangian least-squares isosurface representationprovides a spatio-temporal compression of the output of surface techniques that isamenable to parallel reconstruction in real-time and varying levels of detail.Together these contributions lay a foundation of optimal simulation techniquesand playback for high fidelity fluid animation in games and virtual reality. Theyalso provide improvements on the state of the art for more traditional film orientedvisual effects.82Bibliography[1] Corto mesh compression library. https://github.com/cnr-isti-vclab/corto,2017. → pages 78[2] B. Adcock, D. Huybrechs, and J. Martı´n-Vaquero. On the numericalstability of fourier extensions. Foundations of Computational Mathematics,14(4):635–687, 2014. → pages 73[3] C. Anderson and C. Greengard. On vortex methods. SIAM Journal onNumerical Analysis, 22(3):413–440, June 1985. → pages 28[4] R. Angst, N. Thuerey, M. Botsch, and M. Gross. Robust and efficient wavesimulations on deforming meshes. In Computer Graphics Forum,volume 27, pages 1895–1900. Wiley Online Library, 2008. → pages 46[5] K. E. Atkinson. The numerical solution of integral equations of the secondkind. Number 4. Cambridge university press, 1997. → pages 51[6] Atkinson,K. E. . The Numerical Solution of Integral Equations of the SecondKind. Cambridge Univ. Pr., 1997. → pages 17[7] J. Beale and A. Majda. High order accurate vortex methods with explicitvelocity kernels. Journal of Computational Physics, 58(2):188–208, 1985.→ pages 30[8] J. T. Beale and A. Majda. Vortex methods I: Convergence in threedimensions. Mathematics of Computation, 39(159):1–27, July 1982. →pages 28[9] R. Bridson. Fluid Simulation for Computer Graphics. A K Peters, 2008. →pages 26[10] R. Bridson, J. Houriham, and M. Nordenstam. Curl-noise for proceduralfluid flow. ACM Transactions on Graphics (Proc. SIGGRAPH), 26(3):46,2007. doi:http://doi.acm.org/10.1145/1276377.1276435. → pages 2783[11] T. Brochu and R. Bridson. Robust topological operations for dynamicexplicit surfaces. SIAM Journal on Scientific Computing, 31(4):2472–2493,2009. doi:10.1137/080737617. URL http://link.aip.org/link/?SCE/31/2472/1.→ pages 30, 32, 37, 45[12] T. Brochu and R. Bridson. Animating smoke as a surface.Eurographics/ACM SIGGRAPH Symposium on Computer Animation(posters and demos), 2009. → pages 28[13] T. Brochu, C. Batty, and R. Bridson. Matching fluid simulation elements tosurface geometry and topology. ACM Transactions on Graphics (Proc.SIGGRAPH), 29:47:1–47:9, July 2010. ISSN 0730-0301.doi:http://doi.acm.org/10.1145/1778765.1778784. URLhttp://doi.acm.org/10.1145/1778765.1778784. → pages 27[14] T. Brochu, T. Keeler, and R. Bridson. Linear-time smoke animation withvortex sheet meshes. In Eurographics/ACM SIGGRAPH Symposium onComputer Animation, pages 87–95. The Eurographics Association, 2012. →pages iv, 4, 11, 12, 17, 23, 45[15] H. Cheng, L. Greengard, and V. Rokhlin. A fast adaptive multipolealgorithm in three dimensions. Journal of Computational Physics, 155(2):468–498, 1999. → pages 36, 59[16] N. Chentanez and M. Mu¨ller. Real-time simulation of large bodies of waterwith small scale details. In Proceedings of the 2010 ACMSIGGRAPH/Eurographics symposium on computer animation, pages197–206. Eurographics Association, 2010. → pages 46[17] N. Chentanez and M. Mu¨ller. Real-time eulerian water simulation using arestricted tall cell grid. In ACM Transactions on Graphics (TOG),volume 30, page 82. ACM, 2011. → pages 46[18] N. Chiba, K. Muraoka, H. Takahashi, and M. Miura. Two-dimensionalvisual simulation of flames, smoke and the spread of fire. The Journal ofVisualization and Computer Animation, 5(1):37–53, 1994. ISSN 1099-1778.doi:10.1002/vis.4340050104. URLhttp://dx.doi.org/10.1002/vis.4340050104. → pages 28[19] G.-H. Cottet and P. D. Koumoutsakos. Vortex Methods: Theory andpractice. Cambridge University Press, 2000. → pages 11, 28, 3584[20] M. Duff. Quadrature over a pyramid or cube of integrands with a singularityat a vertex. SIAM journal on Numerical Analysis, 19(6):1260–1262, 1982.→ pages 52[21] D. Enright, R. Fedkiw, J. Ferziger, and I. Mitchell. A hybrid particle levelset method for improved interface capturing. Journal of ComputationalPhysics, 183(1):83–116, 2002. → pages 27[22] Exotic Matter. Naiad 0.6. 2011. → pages 27[23] R. Fedkiw, J. Stam, and H. W. Jensen. Visual simulation of smoke. InProceedings of the 28th annual conference on Computer graphics andinteractive techniques, SIGGRAPH ’01, pages 15–22, New York, NY, USA,2001. ACM. ISBN 1-58113-374-X.doi:http://doi.acm.org/10.1145/383259.383260. → pages 27[24] S. Fleishman, D. Cohen-Or, and C. T. Silva. Robust moving least-squaresfitting with sharp features. In ACM transactions on graphics (TOG),volume 24, pages 544–552. ACM, 2005. → pages 69[25] C. Fochesato and F. Dias. A fast method for nonlinear three-dimensionalfree-surface waves. Proceedings of the Royal Society A: Mathematical,Physical and Engineering Science, 462(2073):2715–2735, 2006. → pages44, 45, 51[26] N. Foster and D. Metaxas. Modeling the motion of a hot, turbulent gas. InProceedings of the 24th annual conference on Computer graphics andinteractive techniques, SIGGRAPH ’97, pages 181–188, New York, NY,USA, 1997. ACM Press/Addison-Wesley Publishing Co. ISBN0-89791-896-7. doi:http://doi.acm.org/10.1145/258734.258838. → pages 27[27] A. Fournier and W. T. Reeves. A simple model of ocean waves. InProceedings of the 13th Annual Conference on Computer Graphics andInteractive Techniques, SIGGRAPH ’86, pages 75–84, New York, NY, USA,1986. ACM. ISBN 0-89791-196-2. doi:10.1145/15922.15894. URLhttp://doi.acm.org/10.1145/15922.15894. → pages 46[28] M. N. Gamito, P. F. Lopes, and M. R. Gomes. Two-dimensional simulationof gaseous phenomena using vortex particles. In Proceedings of the 6thEurographics Workshop on Computer Animation and Simulation, pages3–15. Springer-Verlag, 1995. → pages 2885[29] A. Golas, R. Narain, J. Sewall, P. Krajcevski, P. Dubey, and M. Lin.Large-scale fluid simulation using velocity-vorticity domain decomposition.ACM Transactions on Graphics (TOG), 31(6):148, 2012. → pages 45[30] A. Greenbaum, L. Greengard, and G. B. McFadden. Laplace’s equation andthe Dirichlet-Neumann map in multiply connected domains. Journal ofComputational Physics, 105(2):267–348, 1993. → pages 29[31] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations.Journal of Computational Physics, 73:325–348, December 1987. → pages19, 36, 45[32] D. H. Griffel. Applied Functional Analysis. Dover Publications, 1985. →pages 17[33] S. T. Grilli, P. Guyenne, and F. Dias. A fully non-linear model forthree-dimensional overturning waves over an arbitrary bottom. InternationalJournal for Numerical Methods in Fluids, 35(7):829–867, 2001. → pages51, 55[34] C. Groetsch. Integral equations of the first kind, inverse problems andregularization: a crash course. In Journal of Physics: Conference Series,volume 73, page 012001. IOP Publishing, 2007. → pages 51[35] J. L. Hess. Calculation of potential flow about arbitrary three-dimensionallifting bodies. Technical report, DTIC Document, 1972. → pages 34[36] J. P. M. Hultquist. Constructing stream surfaces in steady 3D vector fields.In VIS ’92: Proceedings of the 3rd conference on Visualization ’92, pages171–178, Los Alamitos, CA, USA, 1992. IEEE Computer Society Press.ISBN 0-8186-2896-0. → pages 27[37] G. Irving, E. Guendelman, F. Losasso, and R. Fedkiw. Efficient simulationof large bodies of water by coupling two and three dimensional techniques.ACM Transactions on Graphics (TOG), 25(3):805–811, 2006. → pages 46[38] D. L. James and D. K. Pai. ArtDefo: accurate real time deformable objects.In Proceedings of the 26th annual conference on Computer graphics andinteractive techniques, SIGGRAPH ’99, pages 65–72, New York, NY, USA,1999. ACM Press/Addison-Wesley Publishing Co. ISBN 0-201-48560-5.doi:http://dx.doi.org/10.1145/311535.311542. URLhttp://dx.doi.org/10.1145/311535.311542. → pages 4586[39] T. Ju, F. Losasso, S. Schaefer, and J. Warren. Dual contouring of hermitedata. In ACM Transactions on Graphics (TOG), volume 21, pages 339–346.ACM, 2002. → pages 74, 75[40] Z. Karni and C. Gotsman. Spectral compression of mesh geometry. InProceedings of the 27th annual conference on Computer graphics andinteractive techniques, pages 279–286. ACM Press/Addison-WesleyPublishing Co., 2000. → pages 65[41] T. Keeler. An integral equation method for solving Laplaces equation withRobin boundary conditions. Master’s thesis, Simon Fraser University, 2011.→ pages 17[42] T. Keeler and R. Bridson. Ocean waves animation using boundary integralequations and explicit mesh tracking. In Eurographics/ACM SIGGRAPHSymposium on Computer Animation, pages 11–19. The EurographicsAssociation, 2014. → pages iv, 4, 43[43] O. D. Kellogg. Foundations of Potential Theory. Dover Publications, 1929.→ pages 6, 34[44] B. Kim, Y. Liu, I. Llamas, and J. Rossignac. Advections with significantlyreduced dissipation and diffusion. IEEE Transactions on Visualization andComputer Graphics, 13(1):135–144, 2007. ISSN 1077-2626.doi:http://dx.doi.org/10.1109/TVCG.2007.3. → pages 27[45] D. Kim, O.-y. Song, and H.-S. Ko. Stretching and wiggling liquids. ACMTransactions on Graphics (Proc. SIGGRAPH Asia), page 120, 2009. →pages 28[46] T. Kim, N. Thu¨rey, D. James, and M. Gross. Wavelet turbulence for fluidsimulation. ACM Transactions on Graphics (Proc. SIGGRAPH), pages 1–6,2008. doi:http://doi.acm.org/10.1145/1399504.1360649. → pages 27[47] H. Krishnan, C. Garth, and K. I. Joy. Time and streak surfaces for flowvisualization in large time-varying data sets. Proceedings of IEEEVisualization ’09, Oct. 2009. → pages 27[48] A. T. Layton and M. van de Panne. A numerically efficient and stablealgorithm for animating water waves. The Visual Computer, 18(1):41–53,2002. → pages 4687[49] D. Levin. The approximation power of moving least-squares. Mathematicsof Computation of the American Mathematical Society, 67(224):1517–1531,1998. → pages 68[50] Y. Liu, M. Xue, and D. Yue. Computations of fully nonlinearthree-dimensional wave-wave and wave-body interactions. part 2. nonlinearwaves and forces on a body. Journal of Fluid Mechanics, 438:41–66, 2001.→ pages 44[51] W. E. Lorensen and H. E. Cline. Marching cubes: A high resolution 3dsurface construction algorithm. In ACM siggraph computer graphics,volume 21, pages 163–169. ACM, 1987. → pages 73[52] F. Losasso, F. Gibou, and R. Fedkiw. Simulating water and smoke with anoctree data structure. ACM Transactions on Graphics (Proc. SIGGRAPH),23(3):457–462, 2004. → pages 27[53] A. Maglo, G. Lavoue´, F. Dupont, and C. Hudelot. 3d mesh compression:Survey, comparisons, and emerging trends. ACM Computing Surveys(CSUR), 47(3):44, 2015. → pages 65[54] G. Morton. A Computer Oriented Geodetic Data Base and a New Techniquein File Sequencing. International Business Machines Company, 1966. URLhttp://books.google.ca/books?id=9FFdHAAACAAJ. → pages 21[55] M. Mu¨ller. Fast and robust tracking of fluid surfaces. In Proceedings of the2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation,pages 237–245, New York, NY, USA, 2009. ACM. ISBN978-1-60558-610-6. doi:http://doi.acm.org/10.1145/1599470.1599501. →pages 27[56] M. Mu¨ller, S. Schirm, and S. Duthaler. Screen space meshes. In Proceedingsof the 2007 ACM SIGGRAPH/Eurographics symposium on Computeranimation, pages 9–15. Eurographics Association, 2007. → pages 75[57] R. Narain, J. Sewall, M. Carlson, and M. Lin. Fast animation of turbulenceusing energy transport and procedural synthesis. ACM Trans. Graphics(Proc. SIGGRAPH Asia), 27(5):166:1–166:8, 2008. → pages 27[58] M. B. Nielsen and R. Bridson. Guide shapes for high resolution naturalisticliquid simulation. ACM Trans. Graph., 30(4):83:1–83:8, July 2011. ISSN0730-0301. doi:10.1145/2010324.1964978. URLhttp://doi.acm.org/10.1145/2010324.1964978. → pages 4688[59] J. Park, Y. Seol, F. Cordier, and J. Noh. A smoke visualization model forcapturing surface-like features. Computer Graphics Forum (Proc.Eurographics), 29(8):2352–2362, 2010. ISSN 1467-8659.doi:10.1111/j.1467-8659.2010.01719.x. URLhttp://dx.doi.org/10.1111/j.1467-8659.2010.01719.x. → pages 28[60] S. I. Park and M. J. Kim. Vortex fluid for gaseous phenomena. InProceedings of the 2005 ACM SIGGRAPH/Eurographics Symposium onComputer animation, SCA ’05, pages 261–270, New York, NY, USA, 2005.ACM. ISBN 1-59593-198-8.doi:http://doi.acm.org/10.1145/1073368.1073406. URLhttp://doi.acm.org/10.1145/1073368.1073406. → pages 12, 28, 32, 45[61] D. R. Peachey. Modeling waves and surf. In ACM Siggraph ComputerGraphics, volume 20, pages 65–74. ACM, 1986. → pages 46[62] T. Pfaff, N. Thuerey, and M. Gross. Lagrangian vortex sheets for animatingfluids. ACM Trans. Graph. (Proc. SIGGRAPH), 31(4):112:1–112:8, 2012.→ pages 28, 29, 40[63] F. Ponchio and M. Dellepiane. Fast decompression for web-basedview-dependent 3d rendering, 2015. URLhttp://doi.acm.org/10.1145/2775292.2775308. → pages 78[64] P. G. Saffman. Vortex Dynamics. Cambridge University Press, 1992. →pages 8, 10[65] H. Schechter and R. Bridson. Evolving sub-grid turbulence for smokeanimation. In Proceedings of the 2008 ACM SIGGRAPH/EurographicsSymposium on Computer Animation, 2008. → pages 27[66] A. Selle, N. Rasmussen, and R. Fedkiw. A vortex particle method forsmoke, water and explosions. ACM Transactions on Graphics (Proc.SIGGRAPH), 24(3):910–914, 2005. ISSN 0730-0301.doi:http://doi.acm.org/10.1145/1073204.1073282. → pages 27, 28[67] A. Selle, R. Fedkiw, B. Kim, Y. Liu, and J. Rossignac. An unconditionallystable MacCormack method. SIAM Journal on Scientific Computing, 35(2-3):350–371, 2008. ISSN 0885-7474.doi:http://dx.doi.org/10.1007/s10915-007-9166-4. → pages 27[68] D. Shepard. A two-dimensional interpolation function for irregularly-spaceddata. In Proceedings of the 1968 23rd ACM national conference, pages517–524. ACM, 1968. → pages 6889[69] J. Stam. Stable fluids. In Proceedings of the 26th annual conference onComputer graphics and interactive techniques, SIGGRAPH ’99, pages121–128, New York, NY, USA, 1999. ACM Press/Addison-WesleyPublishing Co. ISBN 0-201-48560-5.doi:http://doi.acm.org/10.1145/311535.311548. → pages 24, 27[70] M. J. Stock, W. J. Dahm, and G. Tryggvason. Impact of a vortex ring on adensity interface using a regularized inviscid vortex sheet method. Journalof Computational Physics, 227(21):9021 – 9043, 2008. ISSN 0021-9991.doi:DOI:10.1016/j.jcp.2008.05.022. URLhttp://www.sciencedirect.com/science/article/B6WHY-4SPC0W8-1/2/fabd09917ff24cd287e22538193d2b6f. Special Issue Celebrating TonyLeonard’s 70th Birthday. → pages 26, 29, 40[71] J. Tessendorf et al. Simulating ocean water. Simulating Nature: Realisticand Interactive Techniques. SIGGRAPH, 2001. → pages x, 3, 44, 46, 52, 60[72] 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, Mar. 1992. ISSN 0196-5204.doi:10.1137/0913035. URL http://dx.doi.org/10.1137/0913035. → pages 54[73] W. von Funck, T. Weinkauf, H. Theisel, and H.-P. Seidel. Smoke surfaces:An interactive flow visualization technique inspired by real-world flowexperiments. IEEE Transactions on Visualization and Computer Graphics,14(6):1396–1403, Nov/Dec 2008. → pages 27[74] S. Weißmann and U. Pinkall. Filament-based smoke with vortex sheddingand variational reconnection. ACM Transactions on Graphics (Proc.SIGGRAPH), 29:115:1–115:12, July 2010. ISSN 0730-0301.doi:http://doi.acm.org/10.1145/1778765.1778852. URLhttp://doi.acm.org/10.1145/1778765.1778852. → pages 11, 12, 28, 33, 45[75] C. Wojtan, N. Thu¨rey, M. Gross, and G. Turk. Deforming meshes that splitand merge. ACM Transactions on Graphics (Proc. SIGGRAPH), 28(3):1–10, 2009. ISSN 0730-0301.doi:http://doi.acm.org/10.1145/1531326.1531382. → pages 27[76] C. Wojtan, N. Thu¨rey, M. Gross, and G. Turk. Physics-inspired topologychanges for thin fluid features. ACM Transactions on Graphics (Proc.SIGGRAPH), 29:50:1–50:8, July 2010. ISSN 0730-0301.doi:http://doi.acm.org/10.1145/1778765.1778787. URLhttp://doi.acm.org/10.1145/1778765.1778787. → pages 2790[77] G. Wyvill, C. McPheeters, and B. Wyvill. Data structure for soft objects.The Visual Computer, August(2), 1986. → pages 73[78] M. Xue, H. Xu, Y. Liu, and D. Yue. Computations of fully nonlinearthree-dimensional wave-wave and wave-body interactions. part 1. dynamicsof steep three-dimensional waves. Journal of Fluid Mechanics, 438:11–39,2001. → pages 44, 51[79] L. Yaeger, C. Upson, and R. Myers. Combining physical and visualsimulation–creation of the planet jupiter for the film “2010”. SIGGRAPHComput. Graph., 20(4):85–93, Aug. 1986. ISSN 0097-8930. → pages 28[80] J. Yu, C. Wojtan, G. Turk, and C. Yap. Explicit mesh surfaces for particlebased fluids. Computer Graphics Forum (Proc. Eurographics), 31(2), 2012.→ pages 27[81] C. Yuksel, D. H. House, and J. Keyser. Wave particles. ACM Transactionson Graphics (TOG), 26(3):99, 2007. → pages 46[82] K. Zhou, Z. Ren, S. Lin, H. Bao, B. Guo, and H.-Y. Shum. Real-time smokerendering using compensated ray marching. ACM Transactions on Graphics(Proc. SIGGRAPH), 27:36:1–36:12, August 2008. ISSN 0730-0301.doi:http://doi.acm.org/10.1145/1360612.1360635. URLhttp://doi.acm.org/10.1145/1360612.1360635. → pages 37, 66[83] Y. Zhu and R. Bridson. Animating sand as a fluid. In ACM Trans. Graph.(Proc. SIGGRAPH), pages 965–972, New York, NY, USA, 2005. ACM.doi:http://doi.acm.org/10.1145/1186822.1073298. → pages 3, 2791
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Surface based fluid animation using integral equations...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Surface based fluid animation using integral equations : simulation and compression Keeler, Todd 2017
pdf
Page Metadata
Item Metadata
Title | Surface based fluid animation using integral equations : simulation and compression |
Creator |
Keeler, Todd |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | This dissertation looks at exploiting the mathematics of vorticity dynamics and potential flow using integral equations to reformulate critical parts of fully dynamic fluid animation methods into surface based problems. These reformulations enable more efficient calculation and data-structures due to the reduction of the simulation domain to the two dimensional fluid surface, rather than its volume. We also introduce a surface compression and real-time playback method for continuous time-dependent iso-surfaces. This compression method further increases the impact of our highly efficient surface-based simulation methods. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-10-02 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
IsShownAt | 10.14288/1.0355876 |
URI | http://hdl.handle.net/2429/63183 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2017-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_november_keeler_todd.pdf [ 5.16MB ]
- Metadata
- JSON: 24-1.0355876.json
- JSON-LD: 24-1.0355876-ld.json
- RDF/XML (Pretty): 24-1.0355876-rdf.xml
- RDF/JSON: 24-1.0355876-rdf.json
- Turtle: 24-1.0355876-turtle.txt
- N-Triples: 24-1.0355876-rdf-ntriples.txt
- Original Record: 24-1.0355876-source.json
- Full Text
- 24-1.0355876-fulltext.txt
- Citation
- 24-1.0355876.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0355876/manifest