@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Science, Faculty of"@en, "Computer Science, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Woo, Alan Wai Lun"@en ; dcterms:issued "2010-01-08T18:57:18Z"@en, "2006"@en ; vivo:relatedDegree "Master of Science - MSc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """Realistic simulation of smoke is used in the special effects industry to produce smoke in both feature films and video games. Traditional simulations utilize uniformly spaced rectangular computational grids to perform the smoke simulation. Various changes had been proposed to improve different aspects of the simulation, including level of details, memory usage and simulation speed. In this thesis, I propose a novel computational grid that improves upon the level of details as well as memory usage. I propose a frustum aligned grid that takes advantage of the viewing camera because details are most important in the area close to the camera. A frustum aligned grid reduces the amount of grid points necessary to cover the whole domain by placing a high concentration of grid points near the camera while having sparse grid points away from the camera. By using a larger number of grid lines in the direction parallel to the camera and fewer grid lines in the direction perpendicular to the camera, high level of details using a smaller amount of memory can be achieved. The grid is logically rectangular and a perspective transformation can map the grid into a spatially rectangular one. These properties enable the use of existing simulation tools with some modifications, thus maintaining the level of speed. Experimental results and comparison with a standard uniform grid demonstrate the practicality and effectiveness of the proposed method."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/17818?expand=metadata"@en ; skos:note "Realistic Smoke Simulation Using A Frustum Aligned Grid by A l a n W a i L u n Woo B . S c , Univers i ty of B r i t i s h Columbia , 2002 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F S C I E N C E i n T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Computer Science) The University of British Columbia A p r i l 2006 © A l a n W a i L u n Woo, 2006 Abstrac t Realis t ic s imulation of smoke is used in the special effects industry to pro-duce smoke in both feature films and video games. Tradi t ional simulations uti l ize uniformly spaced rectangular computat ional grids to perform the smoke simulation. Various changes had been proposed to improve different aspects of the simulation, inc luding level of details, memory usage and simulation speed. In this thesis, I pro-pose a novel computat ional gr id that improves upon the level of details as well as memory usage. I propose a frustum aligned grid that takes advantage of the view-ing camera because details are most important i n the area close to the camera. A frustum aligned grid reduces the amount of gr id points necessary to cover the whole domain by placing a high concentration of gr id points near the camera while having sparse grid points away from the camera. B y using a larger number of gr id lines i n the direction parallel to the camera and fewer grid lines in the direction perpen-dicular to the camera, high level of details using a smaller amount of memory can be achieved. The gr id is logically rectangular and a perspective transformation can map the gr id into a spatially rectangular one. These properties enable the use of ex-ist ing simulation tools w i t h some modifications, thus maintaining the level of speed. Exper imenta l results and comparison w i t h a standard uniform grid demonstrate the pract ical i ty and effectiveness of the proposed method. n Contents Abstract ii Contents iii List of Tables v List of Figures vi Acknowledgements vi i 1 Introduction 1 2 Tradit ional Smoke Simulation 4 2.1 Governing Equations 4 2.2 T h e G r i d 5 2.3 Simulat ion L o o p 5 2.4 Advect ion 8 2.5 Vor t i c i ty Confinement 9 3 F r u s t u m Al igned G r i d M e t h o d 11 3.1 F rus tum Al igned G r i d 11 3.2 Simulat ion Loop 14 i i i 3.3 Divergence 14 3.4 Advec t ion 19 3.5 Pressure solve 20 3.6 Vor t i c i ty 21 4 Exper imenta l results 24 4.1 Source location 25 4.2 G r i d size 27 4.3 Vor t ic i ty 29 4.4 Compar ison wi th t radi t ional s imulat ion 30 5 Conclusion 36 Bibl iography 37 A p p e n d i x A Laplacian Coefficients 38 iv List of Tables 4.1 Parameters for s imulat ion 25 4.2 R u n time of sample simulations w i t h different gr id sizes in seconds per frame and file size of one frame of smoke data in kilobytes. . . . 28 4.3 e and M I N C U R L used i n figure 4.4 30 4.4 R u n time and memory usage comparison of frustum aligned gr id method and t radi t ional method 32 4.5 Compar ison of s imulat ion t ime used i n subroutines between frustum aligned grid method and t radi t ional method using a 64 by 64 by 64 gr id 32 4.6 G r i d spacing of varies gr id sizes 32 4.7 Statistics for figure 4.6 33 v List of Figures 2.1 Computa t iona l domain and staggered gr id alignment 5 3.1 F rus tum parameters 12 3.2 F rus tum aligned grid 13 3.3 A x i s configuration for X?jkY%k 17 4.1 Images obtained from different smoke simulations. Top: a pulse, Second: two sources, T h i r d : moving source, Bo t tom: approaching source 26 4.2 Simulations w i t h different source locations. Top: back half of frus-tum, M i d d l e : middle of frustum, B o t t o m : front half of frustum. . . . 27 4.3 Compar ison between frustum aligned gr id and rectangular gr id . Top: gr id spacing match back of frustum, Bo t tom: grid spacing match front of frustum 29 4.4 Simulations wi th different e and M I N C U R L 31 4.5 Compar ison of simulations wi th different number of gr id lines i n the C. direction 33 4.6 Compar ison of t radi t ional method and frustum aligned grid method. 34 v i Acknowledgements I 'd like to acknowledge my very understanding supervisor Rober t Br idson , for a l l his help and advice. A L A N W A I L U N W O O The University of British Columbia April 2006 v n Chapter 1 Introduction Smoke s imulat ion techniques have made various advancements in recent years to a point where high-quality realistic smoke can be simulated in reasonable time. T h e major breakthroughs come from the introduct ion of the semi-lagrangian method by S t am [8] and vort ic i ty confinement by Fedkiw et al[ l ] . Together, they form the basic backbone of simulation methods that are ut i l ized today. T h e memory requirement is high t radi t ional ly and various techniques have been developed to address this problem. Rasmussen [6] developed a technique to reduce the memory usage by using interpolated two dimensional flow fields instead of a full three dimensional one. T h e two dimensional flow fields are combined wi th a t i led Kolmogorov velocity field to produce large-scale phenomena such as nuclear explosions. T h i s method is inherently l imi ted by its interpolative nature, thus miss-ing smal l scale details that reside i n the physical space between the two dimensional flow fields. T h e octree data structure introduced by Losasso [5] uses large gr id cells i n area of low activities and smaller gr id cells i n areas wi th interesting flow, reduc-ing the overall number of cells used to divide the computat ional domain. However, the definition of interesting flow and cri ter ia for refinement are difficult to define 1 accurately to capture the required details. Recently, Selle [7] combined the vortex particle method wi th an Euler ian gr id to produce highly turbulent effects. T h e vort ici ty form of the Navier-Stokes equations is solved to obtain velocity for advection and vort ici ty particles are used to conserve the vort ic i ty of the flow. A vort ic i ty conserving force is used to drive the gr id vort ic i ty to match the particles. Fe ldman [2] introduces fluid s imulat ion using unstructured tetrahedral gr id and i n combinat ion wi th regular hexahedral gr id . Th i s method allows for simulation i n irregular domains that are difficult w i th t radi t ional grids. We introduce a novel computat ional gr id that addresses the memory issue as well as level of details. A frustum aligned grid is used to take advantage of the v iewing camera. Unl ike octree, our gr id has a defined area of interest, the front of the frustum, where we place high concentration of small gr id cells to capture the details. Progressively larger gr id cells are placed towards the back of the frustum to reduce the number of cells necessary to divide the computat ional domain. However, the structure of the gr id does impose addi t ional overheads in the computat ion, but they can be compromised by cut t ing back the grid dimension perpendicular to the camera. Since the gr id is aligned w i t h the camera, the dimensions parallel to the camera play a more v i t a l role i n mainta ining the level of details and cut t ing back the gr id dimension perpendicular to the camera has min ima l impact on the quali ty of the images. Divergence and gradient operators on the frustum aligned gr id are developed based on the orthogonal decomposition theorems by H y m a n and Shashkov [4]-Chapter 2 w i l l define the t radi t ional simulation method, chapter 3 w i l l out-line the details of the frustum aligned gr id method, chapter 4 w i l l present some 2 experimental results and chapter 5 w i l l conclude and discuss some future work. 3 C h a p t e r 2 Tradi t ional Smoke Simulat ion Smoke simulations have been used to produce realistic looking smoke for many years and they a l l generally follow the framework discussed in this chapter. 2.1 Governing Equations Smoke is composed of t iny soot particles suspended in air, which behaves according to a set of governing equations known as the Navier-Stokes equations. V P f V - r ) Vt + (V • V ) V + — = ^ — ^ + f P P where p denotes density and r denotes viscous stress tensors. Viscosi ty and compressibili ty are negligible in smoke and density can be assumed to be constant; therefore the equation can be simplified into the incom-pressible Euler equations that conserve bo th mass and momentum. Vt + (V- V ) V + V P = / V - v = o where V = (u, v, w) denotes the velocity, P denotes rescaled pressure and f denotes any external forces such as buoyancy. 4 Figure 2.1: Computa t iona l domain and staggered gr id alignment. 2.2 The Grid Tradi t iona l simulations divide the computat ional domain into equal volume cells. Veloci ty and pressure are arranged i n a staggered grid alignment or M A C gr id [3] where velocities lie perpendicular to the faces of a cell and pressures lie at cell centers. See figure 2.1. T h i s arrangement ensures the correct quanti ty w i l l be at the correct location dur ing the computat ion. Average velocity and divergence of velocity w i l l be collocated w i t h pressure, and gradient of pressure w i l l be collocated w i t h components of velocity. 2.3 Simulation Loop Simulat ing smoke requires solving the incompressible Euler equations for each frame of the simulat ion. T h e typica l frame rate of an animation is 24 frames per second. To solve the Eu le r equations, a projection method is used. The equations are broken 5 up into parts and solved i n a few steps. V1 = Vn - {V • V)VAt V* = V1 + fAt V 2 P * = V • V* T/ra+1 _ y* _ y p * T h e first step is to advect the velocity using the semi-lagrangian method discussed in the next section. B y ignoring the pressure term, we can reduce the Euler equations into V* — Vn - _ - + (V-V)V = / and rearrange to get V * = Vn - (V • V ) V A i + fAt which can be broken up into V1 = Vn - (V • V)VAt V* = V1 + fAt A n intermediate velocity field V1 is found through advection and external forces are added to i t . Ex te rna l forces usually include buoyancy and vort ici ty confinement, which w i l l be discussed i n section 2.5. Boundary conditions are then applied to the velocity and pressure. Types of boundary conditions include Dir ichlet , where values are specified, and Neumann, where the normal derivative is specified. In a smoke simulation, typical boundary condi t ion used for pressure is | ^ = 0 for a l l boundaries and typical boundary conditions used for velocities are | ^ = 0 for open boundaries and velocity =• 0 for the ground. 6 After boundary conditions are applied, the correct pressure to force a diver-gence free velocity is computed. yn+l_y* At + V P = 0 yn+1 _ y* + V P A i = 0 T/n+1 + V P A i = V* B y tak ing the divergence of the equation V . Vn+1 + V 2 P A i = V • V* and using the fact that Vn+1 should be divergence free, V • Vn+1 = 0 we can derive the Poisson equation required to solve for the pressure. V 2 P A t = V • V* A t ime scaled pressure P * = P A t is usually used to simplify the computat ion. T h e divergence of a cell can be computed as / 1 \\ / \\ dx d_ dy \\ dz J U V = Ui+ljk - Uijk + Vij + ik — Vijk + Wijk+l - Wijk and the Laplac ian of pressure can be computed as V 2 P * = 6POA. - R. ijk fi—ljk Pi+ljk fij—lk -fij+lk fijk—1 ^ijk+l T h e pressure is then defined by _ V • V*jk + Pi-ijk + Pi+ljk + Pij-lk + Pij+lk + Pijk-l + Pijk+1 \"ijk — 7. Pi Pi ik. — 1 Pi.'. 7 T h e Poisson equation can be solved using methods such as Gauss-Seidel relaxation, mul t ig r id or conjugate gradient. T h e final step in the simulat ion is to update the velocity by subtracting the pressure gradient. yn+l = y* _ V P A i 2.4 Advection Advec t ion of quantities such as velocity and smoke density are carried out using the semi-lagrangian method introduced by Stam. Semi-lagrangian method is a blend between Eu le r i an methods, which keep track of quantities at fixed gr id points, and Lagrangian methods, which keep track of the location of the quantities over t ime. T h i s method aims to eliminate the instabi l i ty and time step restriction in the com-putat ion. In explicit Eule r ian methods, the t ime step is restricted by the spacing of the gr id nodes. Lagrangian methods have no time step restriction but the lack of a regularly sampled gr id poses other challenges. B y using fixed gr id nodes as i n Eule-r ian methods and tracing the quantity backward similar to Lagrangian methods, the semi-lagrangian method is able to combine the best of both classes of methods into a stable and efficient package. Some disadvantages of the method include numerical dissipation from interpolation and failure to conserve mass and energy. To advect a quantity, the velocity field at the specific location has to be computed. Since smoke densities are located at the cell centers, averaged cell center velocities need to be derived. Vi ijk \\ ( \"i.jfc+»i+1.7fc \\ 2 2 Wj.7fc+Wi.7fc+1 This velocity is used to trace back to where the quantity is coming from. X* = X{jk — VijkAt Interpolation of the eight neighboring values is used to update the current value of the quantity. Trilinear or cubic spline interpolations are normally used but any higher degree interpolation can be used to increase the accuracy of the simulation. Advection of velocities is performed in a similar fashion with averaged veloc-ity field defined at the appropriate cell faces. For example, the velocity field at the location of L 7 i + l j f c is defined as / .. \\ i = 2 ^ Ui+ljk ^ V w i + h k + \\ j Ui+ljk Viik+vi-i+lk+Vi+ljk+Vj+lj+lk 4 'U>iik+1>>iik+\\+'U>i+\\ik+'U>i+lik+l 4 2.5 Vorticity Confinement Vorticity confinement is introduced by Fedkiw et al. in 2000 to attack the problem of small scale details. Since the computation involves a fair amount of averaging or interpolating of quantities, small scale details are lost in the numerical dissipa-tion. Vorticity confinement is an effort to reintroduce the lost details by applying a confinement force to areas of high vorticity. The force is proportional to the grid spacing and disappears as the spacing decreases due to increased grid size. This property is consistent with the Euler equations. The first step is to compute the vorticity of each cell using cell center velocity averaged from the faces. V%jk ^ Uijk Vijk \\ W i j k J / ^ij+lk—Wij-lk —Vjjk+1 +Vj jk-1 \\ 2 LO = S/ XV = V Uj j k+1 — Ujj k -1 ~ Wj +1 j k + Wi -1 j k 2 Vi+ljk—Vj-ljk — Uij + lk+Uij-lk 2 J T h e n vort ic i ty location vectors that point toward higher vort ici ty concentration are computed. n = V H IVIwI Final ly , a confinement force is computed as / = eh(Q x u) where h is the gr id spacing. T h e force is then applied to the velocities. + = At ( fjjk+fi-l.jk \\ 2 fV ,fV 2 fjjk+fjjk-1 2 V wijk j \\ T h e parameter e controls the amount of confinement force and it can have significant influence on the resulting smoke appearance. A s e increases the amount of ac t iv i ty i n the smoke increases and the smoke appears to be more alive. B u t i f e is set too high, it w i l l dominate the other quantities and causes the smoke to appear unrealistic or even uncharacteristic. 10 C h a p t e r 3 Frus tum Al igned G r i d M e t h o d Like vor t ic i ty confinement discussed in the last chapter, many other methods exist that t ry to improve the realism or efficiency of a smoke simulation. Examples of such methods include the vortex particle method, octree data structure and unstructured grids, among others. T h e Frus tum Al igned G r i d method is a method that takes advantage of the viewing camera and tries to improve on both realism and efficiency. 3.1 Frustum Aligned Grid Since the viewing camera can only see so much at a part icular instance, we can concentrate our s imulat ion wi th in the area of interest: the frustum. T h e frustum can be described using four parameters: field of view, aspect ratio, near c l ipping plane and far c l ipping plane. F i e l d of view defines the angle to the left and right of vertical that the camera is capable of capturing. T y p i c a l cameras have a field of view of 30 degrees or less. Aspect ratio defines the ratio of w i d t h to height. C o m m o n aspect ratios are 4:3 for television and 16:9 for widescreen display. Near and far c l ipping planes defines the closest and farthest distance the camera could 11 Nearplnne Field of View Figure 3.1: F rus tum parameters see. See figure 3.1. In computer graphics, a perspective transformation is used to transform a frustum into a rectangular domain. Generally, the rectangular domain w i l l have dimension of -1 to 1 in a l l directions. T h e perspective transformation can be defined as perspective — V tan fov 0 0 0 —aspect tan fov 0 0 0 0 0 0 — (farplane+nearplane) 2(farplane)(nearplane) farplane—nearplane farplane—nearplane and its inverse can be defined as - l perspective — tan fov 0 0 0 0 — tan fov aspect 0 0 0 0 0 0 0 farplane—nearplane — (farplane+nearplane) 2(farplane){nearplane) 2(f arplane)(nearplane) / These two matrices are used to convert posit ion between gr id space and real space. A s in t radi t ional simulations, the computat ional domain is d ivided into cells. 12 Figure 3.2: F rus tum aligned grid We divide the domain using gr id lines that w i l l transform to equally spaced grid lines i n the rectangular gr id space. T h e gr id spacing i n a l l directions w i l l increase gradually from the front to the back of the domain. Th i s produces cells w i t h smaller volume near the camera and cells w i t h larger volume away from the camera. Th i s also produces denser gr id nodes near the camera and sparser gr id nodes away from the camera. Velocities and pressure are s t i l l arranged i n a staggered grid alignment but velocities lay along the gr id lines and pressure resides at the gr id nodes. See figure 3.2. T h i s can be view as a shifted version of the standard M A C grid . If a cell is placed around the node, the pressure w i l l lay i n the center and velocities at the faces. In t radi t ional method, this shift would not cause any change i n the simulation process. Since the grid lines are not orthogonal, velocities values are the projection of the real world values onto the directions of the grid lines. T h e axes for any gr id nodes are £ = (1,0,0), r\\ — (0,1,0) and ( = (x,y,z), where x, y and z varies depending on the location of the gr id node. T h e relationship between velocities and 13 real world values can be described as and more usefully Vy ( 1 0 0 0 1 0 \\ x y z j 1 0 0 0 1 0 -x -y 1 v„ since we will need these conversions when performing later computations. 3.2 Simulation Loop The structure of the simulation loop of the frustum aligned grid method is identi-cal to traditional methods. Essentially, it is still solving the incompressible Euler equations using the projection method. But since the physical domain and the com-putational grid is different, computations such as advection, vorticity, divergence and pressure solve need to be redefined. The following sections describe each in detail. 3.3 Divergence Divergence is the rate of change of a control volume. In a frustum aligned grid, the control volume is the volume around a grid node and the flux is related to the velocities on the grid lines. The divergence operator using the frustum align grid is 14 defined as V • V = -N'^LV -N-1 d_d_d_ dx dy dz L\\\\ L\\2 Liz L21 L22 L23 \\ L31 L 3 2 L33 J where TV is the volume of the node, T is the transpose and L is a 3 by 3 mat r ix that relates the formal and natural inner product of velocities. To compute L, we must first define the natural and formal inner product of the domain. T h e first step in deriving the natural inner product of two velocities X and Y is to convert them into real world values Xr and Yr using the mat r ix defined i n the section 3.1. v XT. J X,, x„ \\ / y r \\ ( v \\ V -Xu(x)-Xv(y)+Xu and 1 v Y„ V -Yu(x)-Yv(y)+Yw J where £ = (x, y, z) is the uni t vector of the direction of the grid line. We can then take the inner product of the two values to get (Xr • Yr) = ( X u Xv -Xu(x)-Xv(y)+Xw ^ Yu Yv -Yu(x)-Yv(y)+Yw z J — ( l + f ? ) XUYU + ^1 + XVYV + -pXwYw+ ^(XUYV + XVYU) + -^{XUYW + XWYU) + ^(XVYW + XWYV) T h i s defines the natural inner product at a specific gr id node and axis alignment. To compute the natural inner product of a cell, the natural inner product of the eight corners of the cell are weighted by the volume it occupies and then summed. T h e natural inner product of the domain is obtained by summing the natural inner product of a l l the cells i n the domain. 15 T h e formal inner product of the domain is defined as [X • Y ] = YJUY1vYJW(XUYU + XVYV + XWYW) which sums the inner product over the domain. The formal and natural inner product of the domain is related by {XR • YR) = [LXR • YR] A p p l y i n g the mat r ix L and expanding the equation, we get ( ( Lu L\\2 L13 \\ f Xu \\ T 1 Yu \\ [LXr • Y r ] = L21 L22 L23 Xv Yv \\ V L31 L32 L33 ) \\ Xw ) ) \\ Yyj J — T,uEv'Ew(LiiXuYu + L\\2XVYU + L\\3XWYU + L2\\XUYV + L22XVYV+ L23XWYV + L,3\\XUYW + L,32XVYW + L33XWYW) T h e explici t formulas for L can be derived by comparing the natural and for-m a l inner products. To derive L u , we need to calculate the equivalent of T,uT,vT,wLuXuYu i n terms of the natural inner products. We need to look at the eight axis configura-tions of natural inner products that contain the same XUYU t e rm and weight them according to their volume wi th in the cell. For X^kY^k, the eight axis configurations are 16 !<». -1.0) (0, i . i i ! (O'jivQ} (0.1,0) (1.0,0) (-1.0.0) (1.0.0) (-1.0,0) ;(6Wi.;b.)--. (oj (0.-1.0) Figure 3.3: Axis configuration for X^kY^k. i c V (1,0,0) (0,1,0) (Xiji Viji Zij) (1,0,0) (0,1,0) ( Viji zij) Vol1 x2 (1,0,0) (0,-1,0) (%ij, Diji V o l k ^ + ,2.) (1,0,0) (0,-1,0) (—Xij, —ytj, ~Zij) (-1,0,0) (0,1,0) ( x i + l j i Vi+lj i z i + l j ) ™ * - l u 1 w (-1,0,0) (0,1,0) (—Xi+lj, — JJi+lj, —Zi+lj) (-1,0,0) (-1,0,0) (0,-1,0) (0,-1,0) (xi+ij,Vi+ij, Zi+ij) (—Xi+lj,—Vi+lj, —Zi+lj) , i+y) ™ * u 1 w where Volk and VoZjJ. are the volumes of the front and back half of cells with index k, respectively, and Volk is the volume of cells with index k. See figure 3.3. The rightmost column represents the volume weighted coefficient of the corresponding 17 natural inner product . B y combining the eight terms, we get the explici t formula for L \\ \\ . LnXijk - 2 Vol L I 2 H — + 2, 2 A i j f c i+ij ) W e can derive the other terms i n the mat r ix L i n similar fashion. 'Voll ijk + V ^ r J l ^ ^ p ( ^ y f c + Xij-lk) + z*+l. [Xi+ljk + Xi+lj-lk) T yw _ o / -J--J.7 / yu> , yu> \\ , —Xj+ij ( V o l k - i YW _L V o l k y w ~ ^ z 2 \\ V o l k - i i J k - l -r Volk^ijk J + ^ V o Z f c _ ! + V o Z k ^ i + l j f c L21Xijk L^Xijk LZ$Xijk (Voll-l , Z ' ^ i j y < 7 /\" , y a \"\\ , Xij+lVij+l (yu , y u \"\\ ^ V ^ I T + V ^ J ^ i j f c + A i - l j f c J + z | + ' + Xi+lj+lk) fc-i ^ Voi f c 0 I ~^ ' .7 / ^ ° ^ f e - i y i o I \" y-k yw \\ I ~yi .7+i / • \" \" f c - i yw _ i _ v \" ' t y i o Z 1 z 2 . ^ V o J f c - i ^ - i j f c - l Volk^ijkJ ^ J ^ - ^ V o / f c . i ^ - i j + l f c - l V o i f c ^ i i + l f e ijk Vij+i I V\"'l l y u £ 3 1 * * = 2 22. VoZ° Vo/fc Vol' :fc \\^ijk -A-i-ljk) Volk {^ijk+l + ^ i - l j f c + l j L^Xijk - 2 (^-^J (^ F^ t ( X ^ f c + Xij-\\k) + V ^ t ( X i j ' f c + 1 + Xij-lk+l) L^xTjk = ^ xTjk tj To compute the divergence, the first step is to compute the temporary values of LV. I utemp \\ ( L l l U + L 1 2 V + L l 3 W L2\\u + L22V + L23W y L3iu + L32V + L33W j T h e n the transpose of gradient is applied and finally mult ipl ied by the inverse of the nodes volume. , .temp W temp \\ V - V = -1 temp Li-\\jk temp temp temp temp Hjk Vij-lk ~ Vijk Wijk-\\ ' A T T - <^ temp • ijk VNodek \\ AXk AYk A % _ i AZijk J AX, AY and AZ are the gr id spacing i n X , Y and Z respectively. T h e grid spacing i n the X and Y directions are constant w i th in the same k index while the grid spacing i n the Z direction varies. 18 3.4 Advection Advec t ion is performed using the semi-lagrangian method as before, but a few changes need to be made i n order to accommodate the frustum aligned gr id . T o advect the smoke density at the grid nodes, the first step is to convert the gr id locat ion into a physical location using the inverse of the perspective matr ix . T h e next step is to define a velocity field at the gr id nodes by averaging the velocities along the gr id lines. Since the velocities are located i n the midpoint between grid nodes and the grid nodes are not equally spaced in the £ direction, we need to do a weighted average instead. V^k V w i j k j Ujjk+Uj-ljk \\ 2 Vjjk+Vij-ik 2 ™iik^Zijk-\\+Wiik-\\&Ziik . &Ziik+&Zijk-\\ / T h e node velocity is then converted to real world values and used to trace the smoke density backward. The calculated physical location is converted back into a grid lo-cation using the perspective matr ix . Tri l inear interpolation is then carried out using the gr id locat ion i n grid space. Special treatment is needed for the interpolation in the £ direction because the perspective transformation is not affine i n that direc-t ion. T h e grid location is used to locate which cell the smoke is from. T h e physical location of the cell's corners and the calculated physical location of the smoke in the C. direction are used to compute a weighting factor for the interpolation. The gr id location is sufficient to compute the weighting factor for interpolation in the £ and rj directions. Higher degree interpolat ion schemes can be used to obtain more accurate results. Unl ike t radi t ional methods where the velocities are advected individual ly , the averaged node velocity is advected instead. The process is identical to advecting 19 smoke described above. T h e advected gr id lines using the usual method. Uijk V wijk J 3.5 Pressure solve To solve for the pressure, we need to solve the Poisson equation defined as V 2P = V • V T h e divergence is defined i n section 3.3 and we s t i l l need to define the Laplac ian of the pressure. We can rewrite the Laplac ian as the divergence of the gradient and apply the definition of divergence from above. Since the divergence operator is not as compact as on a regular grid, the Laplac ian operator involves 18 of the 26 neighboring gr id nodes. V 2P = V • VP = - A T ^ L V P node velocity is then averaged back to the / Uiik+Uj+ljk \\ 2 Vjjk+Vjj+lk 2 2 20 E x p a n d i n g the equation we get L n L12 Lis L21 L22 L23 \\ L31 L32 L33 J W V P ^ V P , V P Z M_1AL11\\7Px+L12VPy+L13VPz)i_ k^ _ ( L i i V P x + L i z V P y + ^ i a V f z ) , ^ J V V AXk AXk _ (.L21VPx+L22VPy+L23VPz)i:j_lk _ {L2lVPx+L22VPy+L23VPz)i;jk + AYk ( L 3 l V P x + L 3 2 V P y + L 3 3 V F 2 ) . . f c _ 1 (L31X7Px+L32VPy+L33VPz)ijk AZ, ijk + _ M - \\ ( {L11VP^j-xjk-jL^VPx)ijk (L12VPy)i_ljk-(Ll2VPy)ijk (Li3VPz)i-iik-(.L13VPz)i:jk J V ^ AXk AXk A X f c \"I\" ( J L2lVPx)i ,- l fc-(^2lVPx)i 7 fc , (I-22VP»)i.7--lfc-(I-22VP»)ijfc , (J-23 VPz)ij-Ifc ~(I<23 VP^ijfc , ( ^ l V P x ) i j t - l _ ( £ 3 l V P * ) i j f c , (^32VPy)ijfc-l _ (L32VPy)iik (LssVPzhjk-l _ (L33VPz)ijk-} Further expanding the above equation by replacing the pressure gradient w i t h / yp? . . \\ / Pj+ijk-Pjjk \\ V P y f c = ijk V Pji + lk—Pjjk Pjjk + l—Pjjk AZ. ijk which defines the gradient of pressure along the gr id lines, we can derive the coef-ficient of the gr id node and its 18 neighbors. T h e derivation of the coefficients is included i n the appendix. T h e pressure can then be solved using the standard Gauss-Seidel relaxation or conjugate gradient. 3.6 Vorticity To increase the realism of the simulation, a force similar to vor i t ic i ty confinement is used. T h e idea is to mainta in the level of cur l w i th in the system, thus maintaining 21 the smal l scale details that would have otherwise be damped out. T h e first step is to calculate the cur l for each face of the cells. c u r l i j k j j A r e a i k V X Vijk = c u r l i j k curlTjk j 'i+\\ik&Zi+\\jk—Uiik+\\&Xk+-L-Wijk&Zijk+Uiik&Xk Areayu 3 K V i + \\ j k & Y k - U i i + \\ k & X k - V i i k & Y k + u i : j k & X K Areaf J T h e n i f the magnitude of the cur l of any face is less then a predefined min imum, a force is applied to increase the cur l of the face. The cur l is always perpendicular to the face and the faces along the same gr id lines lay on the same plane. T h e two dimensional confinement force can be calculated for each plane using the normalized gradient of cur l as a weighting term. V x V%k : weightijk = V wei9htzijk w e i g n t l j k - l w e i g h t . . k l A Y k + A Y k + 1 . c u r l i i k + i - c u r l i i k - i \\ &-Zijk+&Zij+\\k fijk+i = -e(curl?jk)(weightljk)AYk+1/2 /y+ifc = e(curlfjk)(weightyijk)AZij+lk/2 Pijk = -e(curl?jk)(weightzijk)AYk/2 f?3k = e{cuTl^k){weightyljk)AZijk/2 22 V x V?jk : we i 9 h t i j k = . h+z H wpinht-, - wei9ht^k, eigntljk - lweight.'jkl fr+ijk = -e(curlVjk)(weigh%jk)AZi+ljk/2 f?jk+1 = e(curl-jk){weightljk)AXk+1/2 f»ijk = -e{curl\\jk){weightxijk)AZijk/2 fijk = e(curlvjk)(weightfjk)AXk/2 V x V*k : weightijk { , -r. \\ I curl™ ..-curlw , .. \\ weightfjk » / — \\ curlTj+ik-curlTj-ik J \\ w e i 9 h t V i j k J • i • weightnk wghtijk = \\ W e i 9 h t : 3 k \\ fi+ijk = e{curl%k){weight*jk)AYk/2 fij+ik = -e(curlfjk){weightlk)AXk/2 fijk = e{curiy3k){weight*jk)AYk/2 f?jk = -e(curl?jk)(weightyAXk/2 Both e and the minimum curl can be used to control the amount of force applied and the activities within the smoke. They have similar effect as the e used in vorticity confinement where values too large can have negative overall effects. 23 C h a p t e r 4 Exper imenta l results T h e experimental simulations are implemented using V i s u a l C + + and O p e n G L . T h e y are run on a laptop w i t h 1.5Ghz Pen t ium M processor and 5 1 2 M B of memory. T h e simulat ion requires a few parameters including field of view, aspect ratio, near plane, far plane, number of gr id lines in the x, y, and z direction, number of iterations for the Gauss-Seidel pressure solver, number of seconds to simulate, buoyancy, e for vorticity, r for opacity and m i n i m u m cur l . T y p i c a l values used i n experimental simulations are summarized i n table 4.1. Pressure at boundaries are set to zero and Neumann conditions are used for velocities. Exper iments are conducted to investigate the effects of various parameters have on the simulation, and a comparison w i t h t radi t ional simulation is also carried out. A hypothetical mushroom cloud where the smoke originate from a smal l open-ing in the ground i n the center of the domain is used to perform the testing, but i n practice, there can be more than one source and they can be located anywhere wi th in the domain. T h e images are rendered using a ray marching technique which s imply accu-mulates the luminance along the £ direct ion. T w o light sources are used to i l luminate 24 Parameter Value fov aspect near plane far plane 15° 1.0 5.0 10.0 128 128 64 10 20 5.0 2.0 0.1 5.0 X M A X Y M A X Z M A X iterations seconds buoyancy e T M I N C U R L Table 4.1: Parameters for simulation. the smoke, one from overhead and one from the left hand side. More advanced tech-niques such as photon mapping can be used to create higher quali ty images. Figure 4.1 shows several sample images obtained from different styles of smoke. Some artifacts are visible due to aliasing from the gr id . 4.1 Source location Since the configuration of the frustum aligned grid places emphasis more on the space at the front of the frustum than the back, the location of the source should be an important factor i n the appearance of the resultant smoke. Figure 4.2 compares the simulations of three different source locations. These simulations use the typical values found i n the beginning of the chapter except the value of e is set to zero. T h i s effectively turns vort ic i ty confinement off and allows the smoke to rise under buoyancy alone. The sizes of the source i n a l l s imulation are the same and they are i n the shape of a square. The first sequence is generated wi th a source located 25 26 Figure 4.2: Simulations w i t h different source locations. Top: back half of frustum, Midd le : middle of frustum, Bo t tom: front half of frustum. at the back half of the frustum, the second sequence is generated w i t h a source located at the middle of the frustum and the last is generated wi th a source located at the front half of the frustum. Since the size of the source is identical, as it moves further away from the camera, it w i l l appear smaller. T h e level of details w i l l also be lowered as less gr id nodes are located at the back of the frustum, but coupled wi th the reduction in the apparent size, the result ing image st i l l maintains a comparable level of details. 4.2 G r i d size T h e grid size is a crucial factor i n the smoke simulat ion. To increase the realism of a s imulat ion, the t r iv ia l solution is to increase the resolution of the grid since the smaller the gr id spacing, the better the abi l i ty to capture the details. The drawback 27 G r i d Size R u n Time(s/frame) F i l e Size(kb) 32 by 32 by 32 0.461 517 64 by 64 by 32 1.342 2057 64 by 64 by 64 2.874 4113 128 by 128 by 32 5.448 8209 128 by 128 by 128 20.920 32833 256 by 256 by 32 41.309 32801 Table 4.2: R u n time of sample simulations w i t h different gr id sizes in seconds per frame and file size of one frame of smoke data in kilobytes. of such approach includes the increase i n s imulat ion t ime and memory requirement. Typica l ly , the s imulat ion t ime varies direct ly w i t h the grid size; doubl ing the grid dimension in a l l directions transforms to eight times the original s imulation time. Table 4.2 lists the s imulat ion t ime for several sample gr id sizes as well as the file size for one frame of smoke data. T h e frustum aligned gr id method is aimed to increase the realism while addressing the t ime and memory issue. Each grid lines i n the direction perpendicular to the camera is projected into a single point at the front of the camera. T h e dimension parallel to the camera can be viewed as the resolution of the image where each pixel is a projected grid line. Us ing this arrangement, a l l gr id nodes are located at a relevant posit ion wi th respect to the camera and no grid nodes are wasted or overlapping. In t radi t ional rectangular computat ional grids, when rendering the resulting image, usually many different gr id nodes w i l l be mapped to the same pixel , thus wasting the computat ional t ime and memory space for these redundant values. T h e gr id spacing i n a frustum aligned grid is not uniform. The spacing at the back of the frustum can be couple times larger than the ones at the front. Yet , the details are not lost from the sparse spacing because the frustum aligned grid accounted for what could and couldn' t be seen by the camera. To cover the same 28 Figure 4.3: Compar ison between frustum aligned gr id and rectangular grid. Top: gr id spacing match back of frustum, Bo t tom: grid spacing match front of frustum. volume occupied by the frustum using a t radi t ional grid, considerably more gr id nodes are needed since i n order to obtain the same level of details, the gr id spacing of the whole domain should match the gr id spacing at the front and extra gr id nodes are needed outside the frustum to get a rectangular gr id . See figure 4.3. To cover the 64 by 64 by 64 grid of the frustum described by the parameters at the beginning of the chapter, a 128 by 128 by 128 uniform gr id w i l l be needed. The amount of ext ra gr id nodes required is main ly depended upon the near and far c l ipping planes of the frustum. 4.3 Vorticity T h e realism of the smoke simulat ion is enhanced by the introduction of vort ici ty force. T h e force is controlled using two parameters: e and M I N C U R L . T h e M I N -C U R L controls the m i n i m u m cur l desirable. A n y face wi th a cur l smaller than 29 e M I N C U R L 1 2 3 4 1 2 5 5 5 5 10 10 Table 4.3: e and M I N C U R L used in figure 4.4 M I N C U R L w i l l have a force applied to it and that force is proport ional to the cur l itself, e controls the amount of force that w i l l be applied. Figure 4.4 shows several simulations w i t h different setting of e and M I N C U R L and table 4.3 lists the values used i n the simulations. T h e typical values at the beginning of the chapter are used for the other parameters. T h e ideal value of e and M I N C U R L depends upon the scene and varies greatly. T r i a l and error may be required to determine the ideal value but in most case, the ideal value is not required to produce a realistic and interesting simulat ion. 4.4 Comparison with traditional simulation A similar implementat ion to the frustum aligned gr id method of the t radi t ional gr id method is used for the comparison. Table 4.4 compares the run t ime and memory usage of the two methods for several different gr id sizes and table 4.5 compares the simulat ion t ime used i n the subroutines of the two methods. A s seen from the data that the frustum aligned grid method is on average half as fast as the t radi t ional method. T h e extra computat ion t ime is used main ly i n the divergence and pressure solver as expected. To br ing the frustum aligned method on par w i th the t radi t ional method in term of computat ion time, we can 30 Figure 4.4: Simulations w i t h different e and M I N C U R L . 31 Frus tum Tradi t ional G r i d Size R u n T i m e M e m o r y R u n T i m e Memory 32 by 32 by 32 0.461 4124 0.260 3972 64 by 64 by 32 1.342 13812 0.982 12132 64 by 64 by 64 2.874 26500 1.903 23456 128 by 128 by 32 5.448 52104 3.876 46120 Table 4.4: R u n time and memory usage comparison of frustum aligned gr id method and t radi t ional method. Subroutine F rus tum Tradi t ional Ini tal izat ion 0.160 0.010 Advec t ion 0.280 1.472 Vor t i c i ty 0.170 0.160 Divergence 0.160 0.010 Pressure Solve 1.502 0.100 Gradient 0.010 0.010 Table 4.5: Compar ison of simulation t ime used i n subroutines between frustum aligned grid method and t radi t ional method using a 64 by 64 by 64 gr id . reduce the number of gr id lines i n the direction perpendicular to the camera by a factor of two. Since the grid lines are concentrated at the front of the frustum and they are not equally spaced, decreasing the number of gr id lines w i l l have a smaller effect on the sampling at the front of the frustum where details are most important . Table 4.6 shows the gr id spacing of the first few grid lines of varies gr id sizes and figure 4.5 compares the resulting images from the different simulations. T h e result indicates that the reduction of gr id lines produces comparable G r i d Size Szi 5Z2 5zz 5z4 5z5 128 by 128 by 32 128 by 128 by 64 128 by 128 by 128 0.0820 0.0400 0.0198 0.0847 0.0406 0.0199 0.0875 0.0414 0.0201 0.0906 0.0419 0.0202 0.0938 0.0427 0.0204 0.3125 0.1563 0.0781 Table 4.6: G r i d spacing of varies gr id sizes. 32 Figure 4.5: Compar ison of simulations w i th different number of gr id lines i n the ( direction. Figure M e t h o d G r i d Spacing G r i d Size Simulat ion T i m e a Tradi t iona l 0.0419 128 by 128 by 128 14.221 b F rus tum 0.0425 64 by 64 by 64 2.874 c F rus tum 0.0211 128 by 128 by 64 9.684 d F rus tum 0.0211 128 by 128 by 128 20.920 Table 4.7: Statistics for figure 4.6. images that suffer slight detail lost and it also leads to lower memory requirement. F igure 4.6 compares the images obtained from tradi t ional methods and frus-t u m aligned grid method. Three different cr i ter ia are used for comparison including equal grid spacing, equal gr id size and equal computat ion time. Table 4.7 lists the statistics for the simulations. W i t h equal gr id spacing at the front of the frustum, the t radi t ional s imulat ion produces better result than frustum aligned method but the simulation t ime and 33 Figure 4.6: Compar ison of t radi t ional method and frustum aligned grid method. 34 memory usage are bo th five times more. Increasing the grid size by a factor of two i n the x and y direction produces comparable result. T h e memory usage is half of the tradit ional 's and the s imulat ion t ime is about 30% less. W i t h equal gr id size, the frustum aligned method produces superior s imulat ion than t radi t ional method at a cost of increased simulat ion time. These results confirm the performance of the frustum aligned grid method i n terms of memory usage and level of details. 35 Chapter 5 Conclusion We presented a frustum aligned grid for realistic smoke simulation that benefits from knowledge of the viewing camera. The novel computation grid addressed the problem of memory usage and level of details by utilizing non uniform frustum aligned grid cells. Significantly less grid nodes are required to cover the compu-tation domain, thus reducing the memory usage and computation time. The grid lines in the direction perpendicular to the camera can be reduced without suffering detail lost. The reduction further reduces the memory usage and they can trans-late to increased grid lines in other directions, thus increasing the level of details. Comparison with traditional simulation proved the validity and effectiveness of the method. The currently fixed frustum is a limitation and a topic for future work. The current configuration can realistically simulate a large scale phenomenon at a fixed viewpoint but a fly-by scenario with a moving camera will require additional modification, possibly in the form of particles. 36 Bibl iography [1] R . Fedkiw, J . S tam, and H . W . Jensen. V i s u a l s imulation of smoke. In Proc. of SIGGRAPH 2001, pages 15-22, 2001. [2] B . E . Feldman, J . F . O B r i e n , and B . M . Kl ingner . A n i m a t i n g gases wi th hybr id meshes. In Proc. of SIGGRAPH 2005, 2005. [3] F . Har low and J . Welch. Numer ica l calculation of time-dependent viscous in -compressible flow of fluid w i t h a free surface. In The Physics of Fluids 8, pages 2182-2189, 1965. [4] J . M . H y m a n and M . Shashkov. T h e orthogonal decomposition theorems for mimetic finite difference methods. SIAM Journal on Numberical Analysis, 36(3):788-818, 1999. [5] F . Losasso, F . G i b o u , and R . Fedkiw. Simulat ing water and smoke wi th an octree data structure. In Proc. of SIGGRAPH 2004, pages ??-?? , 2004. [6] N . Rasmussen, D . Nguyen, W . Geiger, and R . Fedkiw. Smoke simulation for large scale phenomena. In Proc. of SIGGRAPH 2003, pages 703-707, 2003. [7] A . Selle, N . Rasmussen, and R . Fedkiw. A vortex particle method for smoke, water and explosions. In Proc. of SIGGRAPH 2005, pages 910-914, 2005. [8] J . S tam. Stable fluid. In Proc. of SIGGRAPH 1999, pages 121-128, 1999. 37 Appendix A Laplacian Coefficients = (2( £ f e + ^ ) ( 2 + f + fe)(Pi? AA- 1 7\")-^i2vpy_ l j f c -^i2vpy- f c AXFC i rVoll-i i Y£H\\ixi-yyi-^\\p. A X f c A Y f c V y o / f c . ! \"I\" V o i f c A 2 2 _ ^ j n - l j - l f e 1 (Voll-i , ^ W ^ + l ^ i + l ^ p I 1 f Volk-\\ , V o ^ w Z j + i j y j + i j N p 38 ^ 3 V ^ L l , 7 - f c - ^ 3 V ^ - f c 2 ( - X i - i j V M - i V-oiS 2 c - i j - i A / M - n p AXkAZijk^ V i i V o l f c - i J-n-ljfc-1 i 2 / - i i - i 1 - u y \" ' M p AXkAZijk ^~zrT7^Volk ^ i - l j ' f e + l 2 _ / Z * ! ± l i W V o i f c - i V o % ^ P ^ .. + AxkAZiik_i (\"^7ri)(y02fc_i )- pi+iifc-i ' AXkAZijk\\ zf+1. )\\Volk - 1 (Voli-i , K w % - i i i . j - n p XJ — 1 zj — 1 ^AXkAYkWol^^ Volk>y z 2 . + i J - H - l j + l f c ^ 2 2 v / ^ . _ l t - L 2 2 v p * ; . t A y ^ VoZ f c _i + Vol° Volk 2 r + V< A y ^ Volk + A y 2 ( + Vol° Volk ) ( 2 + ^ + f l±l )p ij ) ( 4 + ^ + | t i - r ^ ^ 3 V P - _ l f c - ^ 2 3 V P / , f c A y f c A Y k l ^ I 7 n ^ k - i A Z i j t - i VolkAZijk)^3-lk 2 f-Vij-i \\ ( V o l k - i \\ p AYkAZijk^ ^l£T\"A V o / f c _ i ) r i j - l k - l ^ A y f c A Z i , - f c V z 2 . . , A v o i k i r i J - l ) ! + l -yii-^s ,Vol° , A y f c V ^ T / A y o / ^ i A Z y f c . ! Voi fcA^-fc^U+l fc ,Voll_ 2 / - V i j + i u ^ ' t - n n \" A y f c A z i 3 f c _ x v - i j r ^ - n vrffc-! ^i j+ifc- i , v o i ° . _2 / - y i j + i ' A y f c A Z i J f c 39 2 / - Z i j w V o i f c - i \\ p j 2 / -Xjj x / V o ^ \\ p ^ A X f c + 1 A 2 i 3 fc I \"IF A y 0/ f c J-^i-lj fc+1 + A X ^ zl JWolk^AZijk..-, VolkAZijk>^+W ,Voll \" A J t + i A Z l j k ( ^ ) ( v^)Pi+ljfc+l &Zijk-\\ &Zijk __2_(zm)( Y<=i v < )P 1 t A y ^ ^ A V o l ^ A ^ . ! l / o / f c A ^ - f c ^ U - l f c 2 { - y i j \\ / V o l k - i \\ p A y i f c . j A z ^ f c . ^ ^ A i / o i f c - i ^ y - i f c - i + A y f c + 1 A z i 3 : k ( z ^ L ) ( v o t ' > P ^ - l k + 1 i 2 r - y n r V o ^ - 1 i P , ~ r A Y k \\ ~ z ^ ~ ) W o l k . 1 A Z i j k ^ 1 VolkAZijkl^3+^k + A Y k . 1 A Z i j k . 1 V o £ _ i ) P i j + l k - l '-ijk-. - l } -V; . Voll A Y k + 1 A Z i j k AZijk-l &-Z