CGI2015 manuscript No.
(will be inserted by the editor)

Incompressibility-Preserving Deformation for Fluid Flows
Using Vector Potentials
Syuhei Sato · Yoshinori Dobashi · Yonghao Yue · Kei Iwasaki ·Tomoyuki Nishita

Abstract Physically-based fluid simulations usually require expensive computation cost for creating realistic
animations. We present a technique that allows the user
to create various fluid animations from an input fluid
animation sequence, without the need for repeatedly
performing simulations. Our system allows the user todeform the flow field in order to edit the overall fluid
behavior. In order to maintain plausible physical behavior, we ensure the incompressibility to guarantee themass conservation. We use a vector potential for representing the flow fields to realize such incompressibilitypreserving deformations. Our method first computes(time-varying) vector potentials from the input velocity field sequence. Then, the user deforms the vector potential, and the system computes the deformed velocity
field by taking the curl operator on the vector potential.
The incompressibility is thus obtained by construction.
We show various examples to demonstrate the usefulness of our method.
Keywords flow deformation · incompressibility ·vector potential

Syuhei Sato
UEI Research, Tokyo, Japan
E-mail: syuhei.sato@uei.co.jp
Yoshinori Dobashi
Hokkaido University, JST CREST, UEI Research, Sapporo,Japan

1 Introduction
Physically-based simulations are ubiquitous for creating
realistic fluid animations. Their expensive computation
cost is, however, usually a common drawback. Because
a number of simulations with diﬀerent parameters haveto be tested before the desired animations are obtained,
the total simulation time could be too long.
In the production environment, artists often create
various fluid animations by deforming the simulated
flow fields or density fields. For example, one may use
procedural methods to deform the flow fields [8, 14], ormay directly deform the density fields. Although theseapproaches allow for generating various stylized animations without the need for repeatedly running the simulations, the results are typically not guaranteed to satisfy the physical conservation laws; violating the massconservation can easily lead to noticeable visible artifacts such as increasing/decreasing fluid mass.
We present a method for deforming existing fluid
flow fields while preserving incompressibility (i.e., massconservation), by using vector potentials. Our method
first computes a (time-varying) vector potential Ψ fromthe input velocity field sequence. Then, the user deforms the vector potential, and the system computes
the deformed velocity field by taking its curl. Because
we have the identity ∇ · ∇×Ψ = 0 for any vector potential Ψ, our technique guarantees the incompressibility
by construction. By using our method, the user candeform existing flow fields to eﬃciently create various
fluid animations.

Yonghao Yue
Columbia University, NY, USA
Kei Iwasaki
Wakayama University, UEI Research, Wakayama, Japan

2 Related WorkTomoyuki Nishita
UEI Research, Hiroshima Shudo University, Tokyo, Japan

Solving the Navier-Stokes equations. Since Stam
had introduced an unconditionally stable solver [21] for2Syuhei Sato et al.

the Navier-Stokes equations, many methods have been
proposed for simulating various fluid phenomena [2,5–7, 10,15]. Usually, parameters need to be adjusted toobtain a desired animation; a naı̈ve adjustment would
require repeatedly performing simulations with diﬀerent parameter sets, possibly leading to an unacceptableworkload.
Fluid control. Several control methods have been presented for creating fluid animations with the desired
shapes (e.g. characters, logos) [4, 20,24], or for enabling
the fluids to move along user specified 3D curves (paths)[13]. These methods control fluid motions by adding external forces to guide the fluids. Using these methods,
the user needs to typically wait for multiple simulations
to finish, in order to create various fluid animations.
Our method oﬀers an alternative way to create various
fluid animations, without the need for re-running any
simulation.
Model reduction. Model reduction methods [23,25]
prepare many sets of velocity fields obtained by simulating fluids with various parameters and initial conditions. Then, principal component analysis (PCA) isapplied to the velocity fields to generate basis functions.
At the cost of expensive precomputation, the flow field
can be updated eﬃciently by computing the NavierStokes equations in this basis space. Kim et al. [12]presented a method for eﬃcient re-simulation with different parameter settings, by applying PCA to a single
set of velocity fields to generate the basis functions. In
these methods, the user is limited to creating flow fields
that are represented by a linear combination of the basis
functions. In addition, memory consumption for solving
PCA can be a bottleneck for a large database.
Procedural methods. Procedural methods can beused to generate various flows with relatively low computation cost. For example, for editing fire animations,
methods (e.g., [8,14]) have been presented to enable
the user to deform the fire via a curve representing its
route. Pighin et al. [17] proposed a method for editing
simulated flow fields via advected radial basis functions.
Since these procedural methods typically do not take
into account physical constraints like the incompressibility, unrealistic results can be produced.
Flow deformation. By using scalar stream functions,
we have previously shown that 2D flow fields can be
deformed while guaranteeing the incompressibility [18].
In this paper, we show that 3D flow fields can be deformed in an incompressible way as well, by using vector
potentials.

Fig. 1 Overview of our method. For the sake of visual clarity,each 3D vector field is illustrated as a 2D vector field.

Computing vector potentials. The computation ofvector potentials from velocity fields appears in computational fluid dynamics as well. For example, in thevortex methods (e.g., [3, 9, 16, 26, 27]), one needs to convert the vorticity field ω = ∇×v back to the velocity
field v, which can be regarded as the vector potential of
ω. To our knowledge, our method is the first that uses
vector potentials for fluid editing. We discuss boundary
conditions and equations that are appropriate for our
problem.
3 Input and Assumptions
Our input is a set of incompressible velocity fields v(t)(i.e., ∇ · v(t) = 0), where t = 0, 1, · · · , T − 1 represents
the frame count, and T is the number of frames of the
input velocity field. We assume that v(t) is suﬃciently
smooth and is defined in a simply connected closed domain Ω ⊂ R3 (i.e., Ω is genus 0). In addition, we cannaturally assume v(t) is tangential to ∂Ω:(v(t)) · n = 0

at the boundary ∂Ω,(1)where n is normal to the boundary.
In addition, we assume the fluid, represented as adensity field, is advected along v(t). Our method generates various animations by replacing the velocity field,
and then advecting the fluid along the new velocity
field. Henceforth, we omit the notation t for brevity.
4 Overview
Instead of manipulating the velocity or density fields
directly, we deal with the vector potential to guarantee

Incompressibility-Preserving Deformation for Fluid Flows Using Vector Potentials

3

the incompressibility of the deformed flow. In a simple
connected closed domain, the Helmholtz-Hodge decomposition theory [1,22] states that a suﬃciently smooth
vector field v can be decomposed asv = ∇×Ψ + ∇p,(2)

where Ψ is a vector potential, p is a scalar field. ∇×Ψ
is tangential to the boundary of the domain:(∇×Ψ) · n = 0

at the boundary,(3)

and ∇p is perpendicular to the boundary.
Since our input velocity field v is incompressible, wehave ∇p = 0 (as in Appendix A), and hencev = ∇×Ψ.

(4)

Because the identity ∇ · (∇×Ψ) = 0 is satisfied for an
arbitrary vector field Ψ, the curl of the vector potential is always incompressible. Therefore, if we manipulate the fluid via the vector potential, and then convert
the vector potential back to the velocity field, we canguarantee the incompressibility of the deformed flow by
construction.
Fig. 1 shows an overview of our method. In the preprocess, the velocity field v at each time step is converted into the vector potential Ψ. At runtime, the user
deforms the spatial domain to obtain the deformed vector potentials Ψ̃. Then, the curl operator is applied toΨ̃ to generate the deformed velocity field ṽ. Finally, thedensity fields are advected along ṽ.

Fig. 2 Deformation using a control path.

g and h, we have F (g + h) = g̃ + h̃ = F (g) + F (h).
Then, the deformed vector potential will have the form
F (Ψ) = F (∇×Φ + ∇s) = F (∇×Φ) + F (∇s). For an
arbitrary deformation, ∇×F (∇s) is generally non-zero,
because F (∇s) might no longer be a curl free field. We
wish to pin down ∇s such that ∇×F (∇s) is coherent
in time so that it does not cause a fluctuation in the
resulting velocity field. Specifically, we enforce ∇s = 0,
by setting
Ψ·n=0

at the boundary, and

(5)

∇·Ψ=0

for the entire domain.

(6)

To compute the desired vector potential Ψ from thevelocity field v, we first apply ∇× to the both sides ofEq.(4) and obtain∇×(∇×Ψ) = ∇×v.

(7)

Under our setting, the solution to Eq.(7) satisfies Eq.(4)as in Appendix B. Then, from the identity∇×(∇×Ψ) = ∇(∇ · Ψ) − ∇2 Ψ,(8)

we have
5 Computing Vector Potentials∇(∇ · Ψ) − ∇2 Ψ = ∇×v.

Although the input incompressible vector field v canbe related to a vector potential Ψ via v = ∇× Ψ, asin Eq.(4), Ψ is usually not unique. Suppose that Ψsatisfies v = ∇× Ψ, then any vector potential Ψ′ =
Ψ + ∇q also satisfies v = ∇×Ψ′ because ∇×∇q = 0
for an arbitrary scalar field q. We discuss how we pin
down this freedom in a way suitable for deforming atime-varying vector field.
First, we observe that according to Eq.(2), the vector potential itself can be decomposed as Ψ = ∇×
Φ + ∇s, where Φ is a vector field and s is a scalar
field. Now, suppose that Ψ is deformed via a map F ,which, for example, maps each position X to a newposition x. Under this deformation, a vector valued
function g(X) will be transformed to g̃(x) = g(X)1 ,which we write F (g) = g̃. For vector valued functions

Finally, substituting Eq.(6) into Eq.(9), we obtain∂x
1
Or, we could consider g̃(x) = ∂X
g(X) if we want to apply
local rotation, stretch and shear to reorient the vector valued
function.

−∇2 Ψ = ∇×v.

(9)

(10)

We solve this Poisson equation with the boundary conditions Ψ · n = 0 (5), (∇×Ψ) · n = 0 (3) and ∇ · Ψ = 0(6) to obtain the vector potential Ψ. As in Appendix C,
we show that Eq.(10) is equivalent to Eq.(7) under our
setting. Hence the solution to Eq.(10) together with theboundary conditions satisfies Eq.(7) and in turn Eq.(4).
To obtain numerical solution to Eq.(10), we used thebiconjugate gradient stabilized method (BiCGSTAB).

6 Generation of Deformed Velocity Field
In our method, we assume the input velocity fields are
given as a sequence of voxel data. We represent the vector potentials using a grid as well. The deformation isapplied to the grid directly. We can in general use any

Syuhei Sato et al.
deviation in total mass (%)

4

input flow

divergence-free aware

without divergence-free aware

80
60
without divergence-free aware

40
20

divergence-free aware

input flow

0
10
-20

20
30
frame count40

50

Fig. 3 The importance of divergence-free aware deformation. The top left insets show the underlying grid of the domain. The
top right insets show the divergence of the velocity field in the region marked by the red dashed rectangular: blue, green, redcolors indicate negative, zero, positive divergence, respectively. The deviation of the total mass is computed as (mt /m0 − 1),
where the total mass m is computed by integrating the density over the entire simulation domain.

deformation method developed for images or meshes forthis purpose. However, if deformations that cause foldovers are applied, large velocities might be generated
because the fold-overs cause discontinuities in the vector fields. In addition, unintentional drastic changes inthe velocity values might occur if the degrees of deformations are too large, because the deformed grid in thelargely deformed region could be under resolved.
To compute the new deformed vector potential Ψ̃,
defined in the world coordinates, we first resample Ψstored in the deformed grid to the cartesian grid storingΨ̃, and then apply the local rotation corresponding to
the deformation, in order to get the correct orientations
for the vector potentials. Next, the deformed velocity
field ṽ is obtained by applying the curl operator to Ψ̃:ṽ = ∇× Ψ̃. The resulting velocity field always satisfies
the divergence-free condition since ∇ · (∇× Ψ̃) = 0.
7 Results
We used a desktop PC with an Intel Core i7 3930K
CPU, 32GB memory to compute all the examples. To
render 3D results shown in Figs. 5, 6 and the right image
in Fig. 7, we used the physically-based renderer “Mitsuba” [11]. We used control-points or control-paths to
specify the deformation in our system (Fig. 2 shows anexample of using a control-path). The grid is deformed
according to the control-points or control-paths using
a method based on moving least squares [19]. In addition, the deformation is executed on a 2D plane, andthis deformation is then applied to each plane perpendicular to the y direction of the 3D grid (Fig. 2). Videos
corresponding to the following examples can be foundin the supplementary material (Online Resource 1).
The importance of divergence-free aware deformation. The ability to preserve the incompressibilityof the velocity fields is important to conserve the total
mass. We see this importance in Fig. 3. For the sakeof visual clarity, we used a 2D flow field as demonstration. The 2D analog of the 3D vector potential is ascalar function called stream function, computed fromthe 2D velocity field using our previous work [18].
In this example, a non-zero density field is advectedalong a time-varying velocity field generated using a
2D fluid simulation, where constant upward velocities
are set at the center bottom of the simulation space todrive the flow. No density is added or removed at runtime. If we use the moving least squares to deform thevelocity field directly without using the stream function (denoted as “without divergence-free aware”), theresulting deformed velocity field will have non-zero divergence, which can in turn cause the total mass todeviate largely from its initial value. In contrast, by
working with the potential, it is possible to guarantee
the incompressibility all the time, which significantly
reduces the mass deviation. We observed the same tendency in deforming 3D flow fields.
Comparing the edited results against simulated
results in 3D. Fig. 4 shows a comparison of results
created by our method and fluid simulations. Fig. 4(a)shows the original input smoke animation. Figs. 4(b’)
through (d’) are created by the fluid simulations, where
uniform external forces are applied at the left boundaryof the simulation space to make the smoke flow to theright. Figs. 4(b) through (d) are the deformed results
obtained using our method, showing that we can mimic
such flows by using our deformation, and the results
are fairly close under small deformations. Moreover, we
see that the deformed animation is still plausible for a
relatively large deformation in this setting.
3D fluid editing. Next, we show applications of our
flow deformation to 3D fluid editing. The timing information is summarized in Table 1. The computation
time of our method at runtime is about 50.0 times fasterthan the simulation, while the precomputation is only

Incompressibility-Preserving Deformation for Fluid Flows Using Vector Potentials
(a)

(b)

(c)

input flow

(d)

(b’)

deformed flows

5
(c’)

(d’)

simulated flows

Fig. 4 Comparison of results generated by our method and a fluid simulation.

(a)

(b)
input flow

(c)

(d)

deformed flows

Fig. 5 Various stylized animations (b) to (d) from a single
input flow (a). The insets show the underlying grid of thedomain.

2.0 times longer than the simulation for generating the
input flow. Hence the trial-and-error process for creating various animations with our method is much more
eﬃcient than re-running simulations.
From a single input flow (Fig. 5(a)), we can create
stylized animations, namely “slender” (Fig. 5(b)), “fat”(Fig. 5(c)) or “tall” smokes (Fig. 5(d)), by squashing
the grid horizontally, vertically, or by vertically stretching the grid, respectively.
In Fig. 6, we created a scene with rising smokes frommultiple chimneys, by using the results from Figs. 4(b)through (d). With our method, we can eﬃciently create
synthetic scenes with multiple varied fluid animations
from a single simulated dataset.
Fig. 7 shows a stylized smoke from a magic lamp.
The input flow is a vertically rising smoke (Fig. 7(a)),
generated by using a fluid simulation. The source of thesmoke is located at the center-bottom of the simulation
space. By deforming the flow using our method, we cancreate a swirling smoke motion as in Fig. 7(b), which
can be in turn used to generate the image on the rightin Fig. 7, for an interesting smoke motion.
Limitations. Since we deform the vector potential instead of the velocity field, the resulting deformation in
the velocity field can be non-intuitive. For example, because the x−, y− and z− components of the vector po-

Fig. 6 Examples of smoke animations rising from multiple
chimneys using the results from Figs. 4(b) through (d).
Table 1 Computation times per frame measured in seconds.
Tns is for the fluid simulation. Tp and Tr are for preprocess
and runtime in our method, respectively.
Fig.
4
5, 7

grid resolution
256 × 128 × 384
192 × 192 × 512

Tns
146
234

Tp
298
407

Tr
2.9
4.4

tential are coupled through the expression v = ∇×Ψ,
deforming the vector potential in a plane will usuallyalso result in a change in the velocity field in the direction perpendicular to the plane. In addition, too large
deformation can result in a non-intuitive flow as well.
We wish to quantify the satisfaction of the deformation
with respect to the degree and type of the deformationin the future.

8 Conclusions and Future Work
We have presented a method for deforming fluid flow
fields while preserving the incompressibility. The divergence free condition is satisfied by construction: byworking with the vector potentials and converting themback to the velocity fields, we can automatically obtain
the incompressibility. In addition, we have presented
a method for pinning down the freedom in the vector

6Syuhei Sato et al.

(a) input flow(b) deformed flow

Fig. 7 Example of a smoke animation rising from a magic
lamp.

potentials, such that it does not contain any curl freecomponent.
As future work, we wish to identify an intuitive
range of deformation. In addition, we are planning todevelop a method that can deform fluid flow fields whileconserving the momentum as well.

Acknowledgements
This work was supported in part by the JSPS Postdoctoral Fellowships for Research Abroad.

References
1. Bhatia, H., Norgard, G., Pascucci, V., Bremer, P.: The
Helmholtz-Hodge decomposition–a survey. IEEE Transactions on Visualization and Computer Graphics 19(8),1386–1404 (2013)
2. Bridson, R.: Fluid Simulation for Computer Graphics.
AK Peters (2008)
3. Cottet, G.H., Koumoutsakos, P.D.: Vortex Methods:
Theory and Practice. Cambridge University Press (2000)
4. Fattal, R., Lischinski, D.: Target-driven smoke animation. ACM Transactions on Graphics 23(3), 439–446(2004)
5. Fedkiw, R., Stam, J., Jansen, H.W.: Visual simulation ofsmoke. In: Proceedings of ACM SIGGRAPH 2001, pp.
15–22 (2001)
6. Feldman, B.E., O’Brien, J.F., Arikan, O.: Animating suspended particle explosions. In: Proceedings of ACM SIGGRAPH 2003, pp. 708–715 (2003)

7. Foster, N., Fedkiw, R.: Practical animation of liquids. In:
Proceedings of ACM SIGGRAPH 2001, pp. 23–30 (2001)
8. Fuller, A.R., Krishnan, H., Mahrous, K., Hamann, B.,
Joy, K.I.: Real-time procedural volumetric fire. In: Proceeding of the 2007 symposium on Interactive 3D graphics and games, pp. 175–180 (2007)
9. Gamito, M.N., Lopes, P.F., Gomes, M.R.: Twodimensional simulation of gaseous phenomena using vortex particles. In: In Proceedings of the 6th Eurographics
Workshop on Computer Animation and Simulation, pp.
3–15. Springer-Verlag (1995)
10. Hong, W., House, D.H., Keyser, J.: Adaptive particles forincompressible fluid simulation. The Visual Computer24(7-9), 535–543 (2008)
11. Jakob,
W.:
Mitsubarenderer(2010).
Http://www.mitsuba-renderer.org
12. Kim, T., Delaney, J.: Subspace fluid re-simulation. ACM
Transactions on Graphics 32(4), Article 62 (2013)
13. Kim, Y., Machiraju, R., Thompson, D.: Path-based control of smoke simulations. In: Proceedings of the 2006
ACM SIGGRAPH/Eurographics symposium on Computer animation, pp. 33–42 (2006)
14. Lamorlette, A., Foster, N.: Structural modeling of flames
for a production environment. ACM Transactions on
Graphics 21(3), 729–735 (2002)
15. Nguyen, D.Q., Fedkiw, R., Jensen, H.W.: Physically
based modeling and animation of fire. ACM Transactions on Graphics 21(3), 721–728 (2002)
16. Park, S.I., Kim, M.J.: Vortex fluid for gaseous phenomena.
In: Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pp. 261–270 (2005)
17. Pighin, F., Cohen, J., Shah, M.: Modeling and editingflows using advected radial basis functions. In: Proceedings of the 2004 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp. 223–232 (2004)
18. Sato, S., Dobashi, Y., Iwasaki, K., Yamamoto, T.,
Nishita, T.: Deformation of 2D flow fields using stream
functions. In: Proceedings of SIGGRAPH Asia 2014
Technical Briefs, Article 4 (2014)19. Schaefer, S., McPhail, T., Warren, J.: Image deformation using moving least squares. ACM Transactions on
Graphics 25(3), 533–540 (2006)
20. Shi, L., Yu, Y.: Taming liquids for rapidly changing targets. In: Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pp. 229–236 (2005)
21. Stam, J.: Stable fluids. In: Proceedings of ACM SIGGRAPH 1999, Annual Conference Series, pp. 121–128(1999)
22. Tong, Y., Lombeyda, S., Hirani, A.N., Desbrun, M.: Discrete multiscale vector field decomposition. ACM Transactions on Graphics 22(3), 445–452 (2003)
23. Treuille, A., Lewis, A., Popovic, Z.: Model reduction forreal-time fluids. ACM Transactions on Graphics 25(3),826–834 (2006)
24. Treuille, A., McNamara, A., Popovic, Z., Stam, J.:
Keyframe control of smoke simulations. ACM Transactions on Graphics 22(3), 716–723 (2003)
25. Wicke, M., Stanton, M., Treuille, A.: Modular bases for
fluid dynamics. ACM Transactions on Graphics 28(3),Article 39 (2009)
26. Yaeger, L., Upson, C., Myers, R.: Combining physicaland visual simulation-creation of the planet Jupiter for
the film “2010”. In: Proceedings of ACM SIGGRAPH’86, pp. 85–93 (1986)

Incompressibility-Preserving Deformation for Fluid Flows Using Vector Potentials
27. Zhang, X., Bridson, R.: A PPPM fast summation method
for fluids and beyond. ACM Trans. Graph. 33(6), 206:1–
206:11 (2014). DOI 10.1145/2661229.2661261. URL
http://doi.acm.org/10.1145/2661229.2661261

A ∇p = 0 for an Incompressible Vector Field v
Taking the divergence of both sides of Eq.(2), we have 0 =∇ · v = ∇2 p. Because v and ∇×Ψ are both tangential to theboundary, we have (∇p) · n = 0 at the boundary. Therefore,
the solution to p is p = const., hence ∇p = 0.

B Equivalence Between Eq.(4) and Eq.(7)Eq.(7) can be rewritten as ∇× (∇× Ψ − v) = 0. Let A =∇×Ψ − v. Since ∇ · ∇×Ψ = ∇ · v = 0, we have ∇ · A = 0. In
addition, since (∇×Ψ) · n = 0 and v · n = 0, we have A · n = 0.
From ∇×A = 0, ∇ · A = 0, and A · n = 0, we have A = 0,
which gives Eq.(4).

C Equivalence Between Eq.(7) and Eq.(10)First, we show ∇·Ψ = 0 given Eq.(10) and the corresponding
boundary conditions. From Eq.(10), we have ∇ · (−∇2 Ψ) =∇ · ∇×v = 0. Since Ψ is suﬃciently smooth, ∇ · (−∇2 Ψ) =−∇2 (∇·Ψ). Hence we have ∇2 (∇·Ψ) = 0. Writing G = ∇·Ψ,
we have ∇2 G = 0. Since ∇ · Ψ = 0 at the boundary, we have
G = 0 at the boundary. Therefore, the solution to G is G = 0
for the entire domain, which yields ∇ · Ψ = 0. Next, from∇ · Ψ = 0, Eq.(8) and Eq.(10), we obtain Eq.(7).

7