A Hyperbolic Geometric Flow for Evolving Films and Foams
SADASHIGE ISHIDA, Nikon Corporation and The University of Tokyo
MASAFUMI YAMAMOTO, The University of Tokyo
RYOICHI ANDO, National Institute of Informatics
TOSHIYA HACHISUKA, The University of Tokyo

Fig. 1. Various configurations of soap films produced by our method; falling bubbles on a film (le�), films spanning a twisted wire (center), and a foam (right).
Our geometric formulation robustly handles complex structures of films and foams such as non-manifold surfaces with topology changes.
Simulating the behavior of soap �lms and foams is a challenging task. A direct numerical simulation of �lms and foams via the Navier-Stokes equations
is still computationally too expensive. We propose an alternative formulation inspired by geometric �ow. Our model exploits the fact, according toPlateau’s laws, that the steady state of a �lm is a union of constant mean
curvature surfaces and minimal surfaces. Such surfaces are also well known
as the steady state solutions of certain curvature �ows. We show a link
between the Navier-Stokes equations and a recent variant of mean curvature �ow, called hyperbolic mean curvature �ow, under the assumption of
constant air pressure per enclosed region. Instead of using hyperbolic mean
curvature �ow as is, we propose to replace curvature by the gradient of
the surface area functional. This formulation enables us to robustly handle
non-manifold con�gurations; such junctions connecting multiple �lms are
intractable with the traditional formulation using curvature. We also add
explicit volume preservation to hyperbolic mean curvature �ow, which infact corresponds to the pressure term of the Navier-Stokes equations. Our
method is simple, fast, robust, and consistent with Plateau’s laws, which are
all due to our reformulation of �lm dynamics as a geometric �ow.
CCS Concepts: • Computing methodologies → Animation; Physical simulation; Geometry;
Additional Key Words and Phrases: �uid simulation, soap �lm, geometric�ow, Plateau’s problem
This work is supported by Nikon Corporation, under the collaborative research projectwith the University of Tokyo.
Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor pro�t or commercial advantage and that copies bear this notice and the full citationon the �rst page. Copyrights for components of this work owned by others than theauthor(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or
republish, to post on servers or to redistribute to lists, requires prior speci�c permission
and/or a fee. Request permissions from permissions@acm.org.
© 2017 Copyright held by the owner/author(s). Publication rights licensed to Association
for Computing Machinery.
0730-0301/2017/11-ART199 $15.00
https://doi.org/10.1145/3130800.3130835

ACM Reference format:
Sadashige Ishida, Masafumi Yamamoto, Ryoichi Ando, and Toshiya Hachisuka.
2017. A Hyperbolic Geometric Flow for Evolving Films and Foams. ACM
Trans. Graph. 36, 6, Article 199 (November 2017), 11 pages.
https://doi.org/10.1145/3130800.3130835

1

INTRODUCTION

Ephemeral yet complex shapes of �lms and foams have fascinatedmany people over the years. The steady-state shape of soap �lms haslong been a subject of study in di�erential geometry [Powers et al.
2002; Sullivan and Morgan 1996; Taylor 1976]. The pioneering study
by Plateau resulted in a set of laws that describe the steady-stateshape called Plateau’s laws. Named after him, Plateau’s problem[Rado 1930] relates the steady-state shape with area-minimizing
surfaces given �xed boundaries. One of the �rst Fields Medals wasawarded to the study [Douglas 1931] of Plateau’s problem, and it remains an active topic of study in mathematics even today [Ambrosio
2015; Harrison 2014; Harrison and Pugh 2015].
In addition to mathematical studies, computer simulation of �lms
and foams has also been a challenging problem. A commonly usednumerical solver for the Navier-Stokes equations can be computationally very expensive since it involves a simulation of multi�uid phenomena with extremely thin geometry deforming withsurface tension [Saye and Sethian 2013]. Eulerian approaches are
particularly not feasible since discretization may overlook such thinfeatures of �lms. Recent works [Durikovic 2001; Zhu et al. 2014]
thus employ Lagrangian approaches and explicitly track surfacesof �lms. Da et al. [2015] formulated �lm dynamics as surface-only
simulation of �lms using vortex sheets. Nonetheless, robust andvisually pleasing simulation of �lms is still challenging.
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

199:2 •S. Ishida et al.

We propose a new formulation of �lm dynamics that uni�es ageometric view of �lms in mathematics into physics simulation of
�lms. Based on the study in mathematics, we formulate �lm dynamics as evolving area-minimizing surfaces. We observe that �lm
dynamics is closely related to a recently studied geometric �owcalled hyperbolic mean curvature �ow [Dexing et al. 2009; He et al.
2009; LeFloch and Smoczyk 2008]. Hyperbolic mean curvature �ow
is a variant of classical mean curvature �ow. Each point on an
evolving surface under hyperbolic mean curvature �ow is accelerated in the direction of its mean curvature normal. We modify this
hyperbolic mean curvature �ow to be volume preserving, and showthat it is, in fact, equivalent to a simpli�ed form of the Navier-Stokes
equations for �lm dynamics under some plausible assumptions.
While it is tempting to adopt existing methods for mean curvature �ow to solve our model, junctions where multiple �lms meet,
namely Plateau borders, become intractable since mean curvature isunde�ned. We thus propose to replace mean curvature by the gradient of surface-area functional which is identical to mean curvaturein the absence of such borders. While the proof of this equivalence
has been known in di�erential geometry, we are the �rst to utilize
this fact to enable robust handling of otherwise intractable con�gurations in geometric �ow. Since our formulation directly minimizes
surface area, the steady-state solutions of our model conform wellwith Plateau’s laws. To summarize, our technical contributions are
as follows:• Introduction of volume preserving hyperbolic mean curvature�ow.
• A variational approach for surface area minimization and amultiregional volume preserving technique, allowing us to treat
multiple �lms in a uni�ed manner.
• Fast, simple to implement, and accurate numerical solver for�lm dynamics which is also consistent with Plateau’s laws.
We provide various experiments to analyze the properties of ourformulation. Figure 1 shows some of our results. Our work enables
an accurate simulation of �lms and foams that outperforms previous
approaches without additional complexity or computational cost.

2 RELATED WORKS
2.1 Animation of Soap Films
For years, visual simulation of �lms and foams has been extensively
studied in computer graphics. For tiny bubbles, researchers haveinvestigated various Eulerian–particle hybrid approaches since they
are small enough to be well represented by a collection of particles.
The seminal work of Hong et al. [2008] employs the particle levelset method and incorporates escaped particles as bubble particles
to simulate underwater bubbles. Busaryev et al. [2012] proposed touse a volume-preserving weighted Voronoi diagram to approximate
the geometry of foams. Kim et al. [2007] successfully simulated
centimeter-scale bubbles and enabled explicit volume control byextending the regional level set method [Zheng et al. 2006]. Unlike
these works, we focus on cases where deformation of individual
bubble surface is visible.
A more natural representation for thin �lms would be explicit
surface meshes. Da et al. [2015] proposed a novel vortex sheet model
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

for surface-only �lm simulation, where a scalar circulation quantity is attached to surface meshes to drive the whole motion. In
such a surface-based simulation, handling merging, and splittingof �lm surfaces is a non-trivial issue by itself. Durikovic [2001]
introduced a numerical technique to handle such phenomena and
Zhu et al. [2014] further improved this tracker, making an e�cient representation of geometrically complex simulation possible.
Given the success of these methods, we also evolve a mesh-basedon surface-only �lm simulation. We, however, propose a di�erent
formulation using geometric �ow.

2.2

Geometry and Flows

Geometry of �lms has been a long-lasting subject of interest among
researchers across a broad range of scienti�c �elds including mathematics [Almgren and Taylor 1976] and engineering [Brew and Lewis
2003]. For instance, mathematicians examined the steady states ofsoap bubbles, which exhibit di�erent visually aesthetic structures
depending on their initial con�gurations [Struwe 2014]. Such structures are driven by the area-minimization due to the surface tension
force [Boys 1958], and they are successfully applied for designingarchitectures [Argyris et al. 1974]. We refer interested readers to abook [Isenberg 1978] for an overview of the dynamics of soap �lms.
Our formulation is inspired by mathematical studies of �lms andfoams, especially those works on Plateau’s problem.
Plateau’s Problem. The shape of a soap �lm converges to a local
minimum of the area functional. The steady state is conjectured tosatisfy the geometric con�gurations described by empirical rules
called Plateau’s laws. Plateau’s problem mathematically formulates
this property as a minimization problem of surface area for given
boundaries, enclosed volumes, and connection of regions. Since
Plateau’s problem is an important subject in many �elds, it has
been studied over 200 years since Lagrange introduced it in 1760.
Hutchings [2002] determined the area minimized con�guration of
two connected bubbles often called "a double bubble", and provedthat it follows Plateau’s laws. Cases for more than two bubbles arestill unsolved.
Plateau’s problem is di�cult to solve even numerically. Researchers
have investigated various methods to handle speci�c scenarios of respective con�gurations. Using mean curvature �ow, Pinkall [1993]
and Polthier [2002] successfully computed discrete minimal surfaces. Variants of a multiple regional level set method [Merriman
et al. 1994; Mohammad and Švadlenka 2002; Švadlenka et al. 2014]
addressed closed regions without open boundaries. There are sometools [Brakke 1992; Pan et al. 2012] to investigate area-minimized
shapes in general settings by heuristically deforming surfaces, however, they do not guarantee convergence to the local minimum.
All of the above methods utilize the mean curvature �ow or its
variations and our method shares the same spirit.
Simply using a curvature �ow, however, does not fully solve Plateau’s
problem in general settings with constraints of enclosed volumes.
We analyze the main di�culties as follows. First, curvature �ows
are applicable only for a single surface since the surface normal
or the curvature is not de�nable for the intersection of surfaces
called non-manifold junctions. Second, volume preserving �ows
can be de�ned only for a single closed surface. Regional level set

A Hyperbolic Geometric Flow for Evolving Films and Foams • 199:3

Initial state

Evolution by A/ x

Intermediate state

Volume correction

Next stateFig. 2. Overview of a time step. We first evolve the films by applying force on the points toward the negative directions of @ A/@x . The volumes of closed
regions may decrease at this point. We next correct the volumes by pushing points toward the normals. The correction amount of each boundary film isrelated with the pressure di�erence through the film.

methods may preserve volumes but cannot handle open surfaces.
We propose an equation and a computational algorithm that naturally handles non-manifold junctions and preserves the enclosed
volumes of multiple regions. Seen as a geometry processing tool,
our method numerically solves Plateau’s problem for an arbitrary
union of closed and open surfaces.
Hyperbolic mean curvature �ow. Recently, hyperbolic variantof mean curvature �ow named hyperbolic mean curvature �ow is
introduced in di�erential geometry [LeFloch and Smoczyk 2008]. It
directly models a vibrating membrane whose acceleration is given
only by mean curvature. We introduce hyperbolic mean curvature
�ow to computer graphics to simulate evolving �lms and foams.
Since original hyperbolic mean curvature �ow cannot be de�ned fornon-manifold junctions and does not preserve the inner volumes ofclosed regions, we extend it to resolve these issues.
Curvature �ows for �uid dynamics. Recently, curvature �ows
are used in computer animation to describe phenomena caused bystrong surface tension. Thürey et al. [2010] used a volume preserving �ow to express complex behaviors of �uid surfaces such as
crown splash. Zhang et al. [2012] employed mean curvature �ow fora surface-only simulation of small water droplets. Evolving a surface toward the negative gradient of the area functional reduces thesurface area, which is equivalent to evolving under mean curvature�ow. Misztal et al. [2014] utilized the property for simulating multiphase �ow of immiscible �uids. These methods assume that �uid
surfaces do not have non-manifold junctions. Our model also accounts for surface tension force via a geometric �ow, and it extends
to non-manifold geometry.

3

OVERVIEWOur key idea is to formulate dynamics of soap �lms using hyperbolic
mean curvature �ow. Hyperbolic mean curvature �ow is de�ned as
d 2x= H (x, t )n(x, t ),(1)
dt 2
where d/dt, x, , H , and n denote total derivative, position, a constant, mean curvature, and the surface normal, respectively. Similar
to mean curvature �ow, the stationary solutions of this hyperbolic
variant are still minimal surfaces. This property is consistent withthe fact that the steady state of a soap �lm without trapped air is aminimal surface. Hyperbolic mean curvature �ow is also known tomodel an elastic membrane moving under the surface tension force
which closely resembles soap �lms [Yau 2000].
Hyperbolic mean curvature �ow, however, has been studied mainly
in di�erential geometry as a geometric �ow and its application toanimation of soap �lms is unclear. Moreover, Equation (1) is not directly applicable to �lm dynamics since the mean curvature cannot
be de�ned on Plateau borders.
In order to handle non-manifold surfaces, we propose to replace
the mean curvature term with the variational derivative of thearea functional with respect to position denoted by @A(x )/@x. This
derivative is known to be equivalent to the mean curvature normal ofsmooth surfaces, and it naturally extends to non-manifold junctions
such as Plateau borders. To account for incompressible air trappedby soap �lms, we include the pressure term pn. Our �nal model
thus becomes
d 2x
=
dt 2

@A(x )
+ pn
@x

(2)

where the constant is equal to twice the surface tension coe�cient
as we show later. The �rst term on the right-hand side of Equation (2)no longer poses complications associated with non-manifold structures since area is measurable regardless of the topology. As weshow later, our model can, in fact, be derived from the Navier-Stokes
equations, under the assumption of constant pressure per bubble
and the atmosphere. We note that this assumption is commonly
adopted for simulation of bubbles[Busaryev et al. 2012; Patkar et al.
2013; Zhu et al. 2014].
In our algorithm, we �rst numerically evaluate the area derivative
to evolve the surface meshes. Next, we assume that pressure is wellapproximated by a constant per trapped region [Patkar et al. 2013;
Zhu et al. 2014]. Under this assumption, the pressure di�erence acts
only between adjacent regions, conveniently allowing us to perform
volume correction by only solving a small and sparse linear systemof which the degrees of freedom is the number of the enclosed
regions. Finally, we apply a symplectic integrator for time evolution.
Figure 2 illustrates the overall steps of our algorithm.

4

GEOMETRIC FLOW OF SOAP FILMS

We now explain our geometric �ow for �lm dynamics. In this paper,a surface is a continuous map U ! R3 where U is a connected
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

199:4 •S. Ishida et al.

ijFig. 3. Illustration of the extended cotangent matrix. We sum up the contribution of all the Voronoi regions of the faces incident to the edges connecting
vertices i and j. This figure shows one of the three Voronoi regions.

subset of R2 . Films consist of manifold and non-manifold geometries.
A manifold point is a point that belongs only to a single surface.
A non-manifold junction is a point that belongs to two or more
surfaces. Especially, a triple-junction is an intersection of threesurfaces, called a Plateau border for soap �lms. A quad-junction
is an intersection of four surfaces. Almost every quad-junction
appearing in soap �lm dynamics is a point where four Plateau
borders meet.
A constant mean curvature (CMC) surface is a surface with thesame mean curvature everywhere. A minimal surface is a special
case of CMC surfaces with zero mean curvature. A minimal surface
locally minimizes its area, that is, any in�nitesimal smooth changeof the surface increases the area. According to Plateau’s laws, thesteady state of a soap �lm is a union of CMC and minimal surfaces.
For example, a double bubble with di�erent pressures consists ofthree CMC surfaces. The steady state of a �lm between two rings
is a minimal surface also known as a catenoid. In general, a �lm
portion with di�erent air pressure on each side becomes a CMC
surface and equal pressure on both sides leads to a minimal surface.
Our model leverages this geometric relationship between soap �lms
and CMC/minimal surfaces for animation of �lms.

4.1

Variational Derivative of the Area Functional

It is well known that mean curvature is unde�ned on discrete surfaces since di�erentiation needs a continuous setting. Discrete differential geometry circumvents this limitation by discretizing mean
curvature normal Hn as di�erentiation of the area functional, basedon the fact that Hn = @A/@x on a smooth surface (see Appendix A).
We found that this equality Hn = @A/@x naturally extends to evennon-manifold junctions, making our geometric �ow de�nable regardless of con�gurations.
Both Hn and @A/@x have the following properties in common:
for an in�nitesimal smooth change of the surface, the change along
its direction gives the maximal area change, and the magnitude indicates how di�erent the current shape is from the area-minimizingcon�guration. By acting force on each point toward the negative
direction of @A/@x, we can locally minimize the surface area. Therefore, if Plateau’s laws as conjectures describing the geometric conditions of the local minima of the area functional are true, ourformulation is automatically consistent with Plateau’s laws.
For a polygon mesh, we discretize @A/@x as follows. It is known
that mean curvature normal is discretized as H N = M 1 LX where
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

Fig. 4. Regions of a double bubble. We correct the volume lossesof the closed domains 1, 2 by moving points of the boundaries@ 12, @ 1Air , @ 2Air toward the normals. We determine the correction amount for each boundary by solving two equations, Vi =Ar ea (i, Air )d i + Ar ea (i, j )(d i d j ) for (i, j ) = (1, 2) and (i, j ) = (2, 1).

matrices H N and X encode mean curvature normals and positionsof the vertices, M is the mass matrix, and L is the cotangent matrix [Desbrun et al. 1999; Jacobson et al. 2013; Meyer et al. 2002].
Note that L is unde�ned for non-manifold vertices. Replacing L by
our extended cotangent matrix L 0 , we can de�ne @A/@x = M 1 L 0X .
Each entry of L 0 is de�ned by,P
1
8
cot( ) for j 2 N (i)
>
f
>
>
>
P 4A or onoi
0<
Li j = >(3)
Lfor j = i
k ,i ik>>>0
otherwise:
where N (i) is the set of neighborhood vertices of i-th vertex and floops over faces incident to the edge connecting x i and x j , andis the opposite angle of the edge in f , and A or onoi is the area ofVoronoi region (Figure 3). It is a natural extension since Li0 j coincides
with Li j if the edge connecting vertices i and j is manifold. Our extended cotangent matrix is sparse similarly to the original de�nition.

4.2

Volume Preservation

Since hyperbolic mean curvature �ow is not a volume-preserving
�ow, volume loss occurs at each time step, which is incompatible
with soap �lm dynamics. We thus introduce the volume preserving
term pn which, as we show later, is equivalent to the pressuredi�erence through the �lm with the assumption of incompressibility
and constant pressure per region.
Instead of explicitly solving for the pressure term, we correcteach enclosed volume by extending Müller’s technique [2009] for
multiple regions. In the original technique, the volume loss V ofan enclosed region is corrected by moving each point on a surface
along its normal direction. The amount of correction per point isgiven as a constant d = V /A where A is the area of the region.
We extend this technique to handle multiple connected regions
by simultaneously preserving all enclosed volumes. In our setting,the amount of correction for each boundary denoted by d relates
to the pressure di�erence of the incident regions. We �rst solve a
linear system Ad = V where d is a vector of coe�cients related
to the inner pressures, and V = (V0init V0 , . . . , VKinit VK ) is avector of volume losses of K enclosed regions. The matrix A is
8>
Area(i, j)
>
>
>
P
<
Ai j = > Area(i, air ) + k Area(i, k )
>
>
>0
:

(regions i and j are adjacent)(j = i)(otherwise)(4)

A Hyperbolic Geometric Flow for Evolving Films and Foams • 199:5
where Area(i, j) is the area of the surface of the boundary domain
between i-th and j-th closed regions, and Area(i, air ) is the area ofthe boundary domain between i-th closed region and the
R external
�uid such as air. Each area is computed as Area(i, j) = @ ds forijeach boundary surface @ i j . Note that A is always invertible.
Figure 4 shows an example of the boundary domains for two connected bubbles. The amount of correction at each point is determinedby d = d (re ionl ) d (re ions ) where re ionl and re ions are tworegions that the point belongs to. For non-manifold junctions, wetake the average of dn over its neighborhoods. For regions of external air, we use d (air ) = 0. Appendix B gives a proof that thisextension corrects all the enclosed volumes.
The volume correction happens after moving vertices by theacceleration d 2x/dt 2 =@A/@x. This operation is analogous tocomputing the pressure term pn with the �rst-order approximation
for time.

ALGORITHM 1: Time integration

4.3

We use a triangle mesh to discretize �lms as in the previous work [Daet al. 2015]. To handle topology changes, we use LosTopos package [Da et al. 2014], which is currently the only surface tracker
that can handle multi-regional surfaces with non-manifold junctions. Since LosTopos requires regions to be closed and any triangle to be incident to two di�erent regions, we use ghost vertices
and triangles that are out of the simulation. For a scene with separate open surfaces, we index 0, 1, 2, . . . , M for open regions andM + 1, M + 2, . . . for closed regions. We re-order the orientation oftriangles such that their normal vectors are always pointing out
towards the smaller number assigned to each region of air.

Connection to Physics

Although Equation (2) appears to lack any physical interpretation, it
can be derived from the inviscid Navier-Stokes equations (Euler
equations) assuming that the �lm is in�nitesimally thin. The Euler
equations with the surface tension force areDu=Dt

H (x )

n+

1

(5)

rp,

where D/Dt, u, , , H , n, and p denote material derivative, thevelocity, �uid density, the surface tension strength, mean curvature,surface normal, and pressure, respectively. The term (x ) is the
Dirac delta function that acts only on the liquid interfaces, essentially encoding the pressure discontinuity.
Integrating both sides of Equation (5) across a thin �lm givesZZ ✓Du
H (x )
1 ◆dV =n + rp dV ,
(6)
Dt
R
R

where R denotes a small domain crossing the interfaces. The lefthand side of the equation simply encodes the temporal change ofthe velocity of the �lm crossing with the domain. Working with the
Dirac function on the right-hand side yields
Z ✓Z
H (x )
1 ◆2 Hn
1n + rp dV =+rpdV .
(7)RR

We have doubled the mean curvature term since the Dirac function
exists twice (e.g., ( 12 ) and ( 12 ) ) where denotes the thicknessof the �lm) on air-liquid and liquid-air interfaces.
Assuming constant pressure for air, meaning that the pressure
gradient only exists within the �lm, the gradient theorem applied
to the second term of the above equation yields
Z
I
p
pback
1
1
1
rpdV =
pndS = front
n=
pn,
(8)
R

@

R

where @ R , pfront , pback , and n denote an in�nitesimal cylinder
across the boundary surface, the pressure at the front side, the pressure at the back side, and the surface normal of the �lm, respectively.
In the limit of in�nitesimal R , the integration domain R approaches a point, allowing us to drop the integral from the left-handside of Equation (6). By using the relation Du/Dt = d 2x/dt 2 , letting= 1, and de�ning = 2 for simplicity, we arrive at Equation (2).

(1) Move constrained vertices.
(2) Set intermediate velocities of free particles.
u0 = u + t (@ A/@x + f ext ).
(3) Evolve free particles to intermediate positions.
x 0 = x + tu 0 .
(4) Update positions by performing volume preservation,x = x 0 + dn.
(5) Update velocities.
u = u 0 + ( d / t )n.
(6) Update the polygon mesh including topology change by LosTopos.
(7) Determine the velocities of newly generated vertices bybarycentric interpolation.

5 IMPLEMENTATION
5.1 Spatial Discretization

5.2

Time Integration

One time step of our time integration scheme is illustrated in Algorithm 1. We take special care to compute the volume correction
vectors dn for non-manifold junctions where normals are unde�ned. For triple-junctions, we average dn of the incident manifold
vertices. For quad-junctions, we average dn of the incident triplejunctions. For junctions more than quad, we handle them in thesame way as quad-junctions. We, however, have never observed
such a junction in our simulation.
Overall, the total computational complexity is O(NT ) + O(N R ) +O(NV ), where NV , NT , and N R are the numbers of vertices, triangles, and regions respectively. Each notation is de�ned as follows:
O(NT ) refers to the cost of computing the volumes and the surface
areas, O(N R ) cost for solving a sparse linear system, and O(NV )for correcting vertices’ positions. We found that O(N R ) scales close
to proportional to the number of non-zero entries in the matrix.
We apply the �rst-order symplectic Euler method. We performthe integration by computing the intermediate velocity u 0 and use
it to integrate the intermediate position x 0 . Once we have computed
x 0 , we perform volume correction.

6

RESULTSWe performed all of our examples on a Mac Pro with 3.5GHz Intel
Xeon E5. We organize our representative examples by topology
change, volume control, scenarios with external forces, computational timings, surface area transition, convergence accuracy to
Plateau’s laws, and a comparison to a real footage.
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

199:6 •S. Ishida et al.

Frame = 20

Frame = 21

Frame = 35

Frame = 44

Frame = 85

Frame = 145

Fig. 5. Foam of six bubbles. The top le� image is a foam at equilibrium.
A�er a bubble bursts at Frame 21, the foam gradually se�les down to thenext equilibrium state while oscillating.

A foam setup of Figure 5 highlights the ability of our method
to reproduce visually convincing soap �lms and foams with highaccuracy. The simulation took 149 msec per timestep on average.
For all the experiments, we employed constant time steps ranging
from 1 to 10 msec depending on the minimal volume of bubbles.
The volume change of regions was less than 0.1% over the courseof simulation. To visualize �lms, we extended the "thindielectric"
BSDF in Mitsuba renderer [Jakob 2010] to achieve the e�ect of thin
�lm interference. Figure 6 shows an example of a �lm strained by
two rings that is stretching and pinches o� as it contacts in thecenter. With the aid of LosTopos surface tracker, our algorithm
naturally handles such topology changes.
Volume Control. Figure 7 shows an example of explicit volume
control due to the in�ation and the de�ation of a bubble. In thisexample, the surfaces of �lms locally bounce for a while after thevolume control and gradually returns to the equilibrium state of aspherical shape. We believe that such an example can be applied toreproduce the other interesting e�ects of injecting or removing air
inside bubbles through a straw. We perform the volume control bysimply changing the volumes of the target bubbles at the beginningof each time step.
External Forces. We point out that external forces such as gravity
forces are straightforward to incorporate in our model by simply
adding them to d 2x/dt 2 with ( @A/@x ). In either case of volume control (Figure 7) and external forces (Figure 9), our operation
does not violate Plateau’s laws. Figure 8 shows an example of our
method coupled with a �uid simulator via external forces. For simplicity, we employed Stable Fluids [Stam 1999], but our method is
not limited to any �uid solver. At the beginning of each time step,
we �rst advance the �uid solver and compute velocity and pressure
on grids. We then obtain the velocity v and the pressure p of the
air at an arbitrary position x through the trilinear interpolation.
When computing the intermediate velocity u 0 , we simply add the
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

Fig. 6. Wireframe view of a strained film a�ached to two rings being
stretched one another. From le� to right: The film splits into two when
the thickness exceeds a critical point. Surface tension forces induce acceleration and they drive the film to split apart and the films converge to twoflat discs.

Fig. 7. Transformation of film surfaces caused by a volume change. Le�: A
double bubble at equilibrium. Center le�: Volume changing by inflating thele� region and deflating the right region. Center right: Each film component
becomes aspherical immediately a�er the volume change for a short while.
Right: The steady state a�er volume change. The boundary membrane ispushed to the larger region since the smaller region has the greater pressure,hence Plateau’s laws still hold.

pressure of the air by
u0 = u + t!
@A+ p(v · n)n .
@x(9)The �lm vertices receiving forces due to the wind are imposed tohave the negative inner product v · n. Figure 8 shows an example
of a bubble blown by the wind spouted out from the bottom.
In this example, we only accounted for one-way coupling fromthe air to �lm. We believe that we can further two-way couple airand �lm by replacing pressure on the grids inside bubbles with ones
computed during the volume correction phase.

6.1

Timings

We analyzed timings through a range of experiments and observedthat our method yields signi�cant performance improvement overprevious methods. Figure 10 illustrates a comparison with themethod of Da et al. [2015] for the same setups. Each timing includes
mesh processing by LosTopos in common. In this experiment, ourresults were from 3 to 21 times faster than Da et al. [2015]. We, however, emphasize that our governing equation and the scheme arequite di�erent from Da et al. [2015]. This comparison thereby serves
as a supplemental. In this experiment, we accelerated the methodof Da et al. [2015] with the fast multipole method (FMM) to improve
the runtime cost down to O (NV ) [Da et al. 2015].
Table 1 outlines timings of our individual simulation step. Currently, the surface evolution and the topological operations via
LosTopos surface tracker dominate the cost of the calculations. However, our method is not tied with any speci�c surface tracker, hence
we view this issue as orthogonal to our main contributions. We noteA Hyperbolic Geometric Flow for Evolving Films and Foams • 199:7

Frame = 0

Frame = 80

Frame = 160

Frame = 240

Fig. 8. A bubble blown up by the wind. The pressure of the air simulated by
a fluid solver a�ects the bubble as an external force.

that "Others" are mostly spent for conversion operations that bridgedi�erent data types of LosTopos engine and Eigen Library [Guennebaud et al. 2010] and thus can be sped up with further low-level
optimization.
Figure 11 compares the transition of our computational complexity of two di�erent scenarios. The top row examines the rise ofcomplexity due to the increase of vertices by increasing the meshresolution. In the bottom row, we increased the number of regionsof cubes in each dimension. The computational cost of our methodscales closely to linear for both cases.

6.2

Analytical Solution and Numerical Convergence

We measured the time transition of surface area for two examples: adouble bubble and a catenoid (See Figure 13 top ), and observed thatthe numerical solution slowly converges to the analytical solution as
expected. Figure 13 shows our results and the counterpart analyticalsurface areas overlaid with striped lines.
The analyticalp surface area of the steady state of a double bubble isAmin = 27 2/3 V /9 for the volume of each representative region
V . We set up a single bubble vertically split at the center with V =2.0612 and Amin = 14.803. The resulting converged area of oursimulation was 14.813.
Our next example is a catenoid, of which the closed-form solution
is given by

Fig. 9. Steady state of a film spanning a cubic frame. Le�: Without external
forces. Right: Under an extremely strong gravity (exaggerated twenty times).

of Brakke [1992] and Da et al. [2015] for reference. Note that discretization of the surface area for each vertex in Equation (10) of Da
et al. [2015] is not uniquely determined at non-manifold junctions.
The errors in T/Q (Da1) column are examples when each of theareas is simply discretized using the incident triangles per region.
Fei [2017] pointed out that discretization as one-third of the total
area of the incident triangles divided by the number of incident regions gives superior results because it eliminates mesh-dependence
and thereby corrects force imbalances arising from asymmetries in
the incident triangulations [Da et al. 2017]. Applying this discretization, results improved to as in T/Q (Da2) column. We observed that
all the triple/quad-junction errors of our examples, including scenes
with external forces, were less than 2.0 and approached zero as themesh resolution increases.

6.4

Comparison with a Real Footage

Figure 12 shows a comparison with a real experiment recorded
using a high-speed camera [Pucci et al. 2015]. In this experiment, an
immersed soap bubble is slowly merging to a static soap �lm at thebottom. We observe that the result of our model resembles the real
footage for a long duration. The subtle di�erence is due to the lossof kinematic energy caused by the way we use the surface tracker
as we discuss later.

(x, , z) = (a cosh( /a) cos u, a cosh( /a) sin u, ) ,where a, b > 0, u 2 [0, 2 ), and 2 [ b, b] denote arbitrary scalarparameters. The analytical solution of the surface area is given
by a 2 (sinh(2b/a) + 2b/a). We set up a cylinder shaped surface
attached to the two solid rings on its edges with parameters set
a = 0.39565, b = 0.28, and Amin = 1.64895, respectively. The
resulting converged area of our simulation was 1.6511. We also
observed that the surface areas of both examples approach theanalytical solutions as we increase the mesh resolution.

6.3

Plateau’s Laws

Plateau’s laws state that three �lm surfaces meet at a triple-junction
at arccos( 1/2) = 120 and four triple-junctions meet at a quadjunction at arccos( 1/3) ⇡ 109.47 [Ball 2009]. These geometricconditions are essential for visually pleasing simulation of �lms. Table 2 provides the list of errors of angles of triple/quad junctions. We
also provide the errors of the same experiments using the method

7

DISCUSSION

Computational complexity. As we noted before, the computational
complexity of our algorithm is O(NT ) + O(NV ) + O(N R ). This ordercan be approximated simply by O (NV ) since NT is linear to NV in
most cases, and the computational cost of volume correction is muchless than the surface evolution by @A/@x even with a large numberof regions. While another formulation by Da et al. [2015] can bemade to achieve the same linear computational complexity using
the fast multipole method, the computational complexity of ourformulation is intrinsically linear without any additional numerical
technique. In our current implementation, the computation timeis dominated by LosTopos. Since our model is independent of thesurface tracker, we expect that a more e�cient surface tracker canpotentially accelerate the overall computation signi�cantly. We,
however, found that LosTopos is currently the best option for our
purpose.
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

199:8 •S. Ishida et al.

TargetSingle cubic bubble
Double bubble
125 latticed bubbles
Film spanning a frame
Catenoid
Six bubbles

Vertices
1003
1792
6702
2679
1346
4882

Faces
2002
3626
16192
5742
2688
9972

Regions
2
3
126
8
2
7

TPF (Ours)
24.70
68.76
160.22431.64769.17
149.04

TPF (Da)
179.72526.833402.61
1323.91221.252570.16

Ratio
7.3
7.7
21.2
3.1
3.2
17.2

Fig. 10. Comparison of timings per frame for various film examples to Da et al. [2015]. Vertices, Faces, and Regions refer to the number of the initial vertices,
faces, and regions, respectively. TPF is the computational time per frame (ms), and Ratio is the ratio of computational timings TPF (Da) / TPF (Ours). The
figure on the right side illustrates the timings of the results in the list. Films A, B, C, and D are a double bubble, a foam of six bubbles, 125 la�ice-configured
bubbles, and a film spanning a frame respectively.
Table 1. Timings per frame by simulation component (ms). The list is ordered
as follows: the total computational time, the evolution by @ A/@x , volume
correction, mesh processing by LosTopos, and remaining operations.
TargetTotal @ A/@x
VC LosTopos Others
Single cubic bubble
125 latticed bubbles
Six bubbles

24.70
160.22149.04

1.01
2.63
4.04

0.01
1.08
1.02

20.98
148.39135.94

2.71
8.12
8.05

Fig. 11. Complexity analysis by di�erent discrete variables. Up: Timings
due to the increase of the number of vertices on surfaces (maximal number
6291458). Bo�om: Timings due to the increase of number of regions (maximal
number 4913). Notice that both draw closely linear curves.

Fig. 13. Top: Initial states of a double bubble and a catenoid and their steadystates. Middle and bo�om: Surface area transition of films. Striped lines
show their corresponding analytical solutions.
Table 2. List of errors of the angles on non-manifold junctions from Plateau’slaws. In the list, T denotes the root mean square (RMS) of the di�erences between arccos( 1/2) and the dihedral angles of the faces incident to Plateau
borders, and Q denotes the RMS of the di�erences between arccos( 1/3)
and the angles of the Plateau borders intersecting at quad-junctions.

Fig. 12. Time sequence of a soap bubble being trapped in a planar film. Up:
Experiment by Pucci et al. [2015]. Bo�om: Our simulation.

ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

TargetT/Q (Ours)

T/Q (Brakke)

T/Q (Da1)

T/Q (Da2)

Double bubble
Film spanning a frame

1.53
0.60 /0.67

1.40
1.80 /2.27

14.67
4.04 /5.77

3.62
1.84 /3.68

A Hyperbolic Geometric Flow for Evolving Films and Foams • 199:9
Alternative Numerical Solvers. To compute intermediate quantities, we can apply an implicit scheme instead of symplectic Euler [Desbrun et al. 1999]. We observed that the implicit scheme
allows us to increase the time step by a factor of ten at most, but it
was not necessarily numerically more stable in our implementation
due to some additional parameters.
Another alternative approach would be to compute the next position and the velocity by computing @A/@x and pn in a strongly
coupled manner without computing intermediate quantities. Since
our model is not analytically di�erentiable, the Newton-Raphson
method is not applicable. Hence, the bisection method is an option
for seeking the value p that preserves the enclosed volumes. We,
however, expect that such a scheme would increase the computational complexity to as much as O(2N R ), and we view this as too
expensive for a large number of bubbles.
Connection to physics. We showed a theoretical connection between our governing equation and the Navies-Stokes equations.
This is valid only for manifold surfaces at this moment. We believethat generalizing quantities such as surface normal for junctions
extends it to non-manifold cases.
Di�erent scenarios. Da et al. [2015] demonstrated scenarios wherea double bubble is pulled apart and bubbles are generated from asweeping wire ring. We note that our method can also handle suchsettings without di�culties. When bubbles are split, we just need toreassign the region labels of the triangles consisting of the boundary
surface. In our implementation, this process is automatically done
by LosTopos. For the latter scenario, we add a region for the newly
generated bubble and compute its initial volume.
Limitations. In comparison to a mesh based simulation of �uid
surfaces that takes thickness of �lms into account [Zhu et al. 2014],
we regard the �lm to be in�nitesimally thin. Our model thus cannot simulate the moment of burst of bubbles or temporal color
change of liquid surface. We expect that incorporation of liquid
simulation on deformable surfaces would enable simulation of suchphenomena using our model.
When a surface tracker generates a new vertex by edge collapse
or split, we interpolate velocities of neighboring vertices by barycentric interpolation. Alternatives for more accurate interpolations arehigher-order approximations such as spline interpolation and Hermite interpolation. We found that, however, a small amount ofenergy loss occurs with any of the higher-order alternatives mentioned above , because a naive interpolation of vectors with di�erent
directions decreases its magnitude, similarly to triangle inequality
for two vectors. We have not found an e�ective solution which
can preserve local energy around a newly generated vertex. Energy
preserving interpolation of a vector �eld on a triangle mesh is aninteresting topic for future work.
Our current implementation is not unconditionally stable. We
speculate that the maximum timestep size is proportional to the maximal edge length and the inverse of the maximal value of @A/@x. We,
however, have not analyzed the CFL condition or relevant condition
for stability since such an analysis depends on several complex factors such as discretization and the stability of a surface tracker. We

also observed that LosTopos becomes unstable if there are signi�cant changes in a triangle mesh. We thus add some damping to avoid
a sudden change for a large time step. A di�erent discretization stepwith an improved surface tracker can potentially resolve this stability issue. This issue is independent of our theoretical formulationof �lm dynamics as a geometric �ow.

8

CONCLUSIONS

We introduced an extension of hyperbolic mean curvature �ow tosimulate dynamics of �lms and foams. Unlike previous works, westart from the observations made by Plateau to construct a geometric �ow which precisely captures the properties of the steady-stateshapes of soap �lms. We generalized a geometric �ow to consistently handle an arbitrary combination of open boundaries, enclosed
volumes, and connectivity of regions, that are crucial for robust simulation of �lms and foams. We showed that the resulting geometric�ow is closely related with the Navier-Stokes equations. The main
advantages of our model are the computational speed and the convergence of the �lms following Plateau’s laws. Our method alsoenables explicit volume control and the straightforward addition of
external force, allowing for new physics phenomena such as �oating
bubbles blown by the wind and then slowly fall due to the gravity.
For Plateau’s problem, our model provides a uni�ed solution as an
area-minimizer without heuristic operations. Moreover, it solves ageneralized version of Plateau’s problem in the presence of external
forces. We believe that it can open up lots of opportunities acrossdi�erent �elds, similarly to the fact that another generalization of
the Plateau’s problem, called Plateau-Kirchho� problem [Giomi and
Mahadevan 2012] led to numerous follow-up works [Giusteri et al.
2017; Perez et al. 2017]. Since our model can also be seen as a general
geometric �ow, future work includes applications to problems ingeometry processing such as �nding CMC surfaces for architecturaldesign.

ACKNOWLEDGMENTS
The authors would like to thank Sigurd Ofstad for discussion andexperiments during the development of the algorithm, Hisanari Otsu
for customization of Mitsuba renderer, Christopher Batty, Fang Da,
and Raymond Yun Fei for the discussion on the previous work, andJamorn Sriwasansak for the advice on the Thai language in the
supplementary video. This work is supported in part by grants from
Nikon Corporation.

A

MEAN CURVATURE AND VARIATIONAL DERIVATIVE
OF THE AREA FUNCTIONAL

A handful of books on Riemannian geometry and mean curvature�ow [Berger and Gostiaux 2012; López 2013; Mantegazza 2011]provide proofs of the equivalence between mean curvature normal
and the variational derivative of the area functional. However, their
proofs are targeted for a Riemannian manifold with any dimension
and not immediately clear how it relates to our setups. To this end,
we provide a sketch of a proof for a special case of surfaces inthree-dimensional Euclidean space. We �rst �x our de�nitions andnotations. The area functional for a surface S : U ! R3 is de�ned
ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

199:10 •S. Ishida et al.

as
A(S ) =

Z

ds =

Z

U

|Su ⇥ S |dud(10)

where = S (U ). Using a continuous map V : U ! R3 , we cande�ne a family of surfaces S = S + V parametrized by a parameter
2 ( E, E) for a small positive number E. For each , we write
L2 ( ) for the space consisting of square-integrable functions on= S (U ). The variational derivative of the area functional on S
is then de�ned by@Ad⌘ arg max@xdkV k=k

=0

A(S )

(11)

R
where kV k is the total speed de�ned by ( |V | 2ds) 1/2 and k is areal constant. To uniquely de�ne @A/@x,
we here set k = kH k withR
a norm in L2 ( ) de�ned by kH k = ( |H | 2ds) 1/2 .

P���� �� @A/@x = Hn. Let S, V , S be as above. We can decompose V into the normal and the tangential part as V = V ?n + V T t.
The �rst variation formula states
Z
Z
d
ds =
V ? Hds
(12)
d =0

where
and ds are the domain of integrationR and the area element for S . Note that the operation F : , h !
(s)h(s)ds is an
inner product of L2 ( ). Therefore we can apply the Cauchy-Schwarz
inequality to the right-hand side of Equation (12), and obtain,
hH, V ? i  kH k kV ? k.

(13)

Equality holds if and only if V ? = C

H H for a constant C H . Therefore, V 0 = C H Hn maximizes the �rst variation among continuous
maps V : U ! R3 with kV k = C for a positive constant C. By taking
C = kH k, we have C H = 1 and obtaindarg max
A(S ) = Hn.
(14)=0kV k= kH k d⇤

B

PROOF OF VOLUME PRESERVATION

We here prove that moving each point by dn with d as in Section 4.2 corrects all the enclosed volumes for �rst-order approximation. This means,Z1(x + dn) · n 0ds 0 ⇡ Vi + Vi
(15)
3 @ 0iwhere @ 0 = S 0 (U ), n 0 = Su0 ⇥S 0 /|Su0 ⇥S 0 |, and ds 0 = |Su0 ⇥S 0 |dud
are the domain of integration, the unit normal, and the area element
for the variational surface S 0 ⌘ S + dn.

P����. We �rst show
R that the left-hand side of Equation (15) isapproximated to Vi + @d ds. We note that non-manifold junci
tions do not contribute to the integral since the set of non-manifold
junctions is of measure-zero with respect to the area element ds.
We have,Su0 ⇥ S 0===

(Su + dnu ) ⇥ (S + dn )Su ⇥ S + dnu ⇥ S + dSu ⇥ n + ( d ) 2nu ⇥ n⇣⌘
1 + d (A + D) + ( d ) 2 (AD BC) Su ⇥ S . (16)

ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.

In the above, we employed Weingarten equations [Kreyszig 1991]
nu = ASu + BS and n = CSu + DS where the coe�cients are
given by
(S · S )(Su · n) (S · S )(Suu · n)*. u+
1.. (Su · S )(Suu · n) (Su · Su )(Su · n) /// .
|Su ⇥ S | 2 . (Su · S )(S · n) (S · S )(Su · n) /D,
, (Su · S )(Su · n) (Su · Su )(S · n) ⇣⌘Hence, it holds ds 0 = 1 + d (A + D) + ( d ) 2 (AD BC) ds andn 0 = n. Using these relations, we haveZ
1
(x + dn) · n 0ds 0(17)
3 @ 0iZ⇣⌘
1(x · n + d ) 1 + d (A + D) + ( d ) 2 (AD BC) ds.
=
3 @ i
A
*.
+
.. B /// =. C /Regarding ( d ) 2 = 0, the above equation is approximated toZZ
1
1x · nds +d (1 + (A + D)x · n) ds
3 @ i
3 @ i
Z
= Vi +
d ds.
(18)
@

i

In the above, we used
the relationsR that A + D = 2H with meanRcurvature H , and @ Hx · nds = @ ds under the assumptioniithat @ i is a closed surface[Kreyszig 1991].
R
We �nally claim @d ds = Vi by the de�nition of d. It
i
holds
Z
Z
XZ
dds =
di ds +
(di d j )ds@@idi Area(i, air ) +=di Aii +

=

@

j,i

i, airXX(did j )Area(i, j)X

d j Ai j .

X

Vk ik = Vi . (20)j,id j Ai j =

j,i

j

=

X
k

k

Vk

X
j

Ai j A jk1 =

(19)

j

Since Ad = V , we have d j = (A 1 V )j =ing this into Equation (19), we obtain
X
X
X
d j Ai j =Ai j
A jk1 Vk
j

i, j

P

1
k A jkVk . Substitut-k

Hence, the enclosed volume of region i is corrected for �rst-orderapproximation.
⇤

REFERENCES
Frederick J Almgren and Jean E Taylor. 1976. The geometry of soap �lms and soap
bubbles. Scienti�c American 235 (1976), 82–93.
Luigi Ambrosio. 2015. Regularity theory for mass-minimizing currents (after AlmgrenDe Lellis-Spadaro). Calculus of Variations and Geometric Measure Theory (2015),
1–23.
John H Argyris, T Angelopoulos, and Bruno Bichat. 1974. A general method for the
shape �nding of lightweight tension structures. Computer Methods in Applied
Mechanics and Engineering 3, 1 (1974), 135–149.
Philip Ball. 2009. Shapes: Nature’s Patterns: A Tapestry in Three Parts (Natures Patterns).
In Shapes: Nature’s Patterns: A Tapestry in Three Parts (Natures Patterns). Oxford
University Press, Oxford,UK, 97–98.
Marcel Berger and Bernard Gostiaux. 2012. Di�erential Geometry: Manifolds, Curves,
and Surfaces: Manifolds, Curves, and Surfaces. Vol. 115. Springer Science & BusinessMedia, New York, NY, USA.

A Hyperbolic Geometric Flow for Evolving Films and Foams •CV Boys. 1958. Soap-Bubbles, their colors and forces which mold them (1890). (1958).
Kenneth A. Brakke. 1992. The Surface Evolver. Experimental Mathematics 1, 2 (1992),
141–165.
JS Brew and WJ Lewis. 2003. Computational form-�nding of tension membrane
structures-Non-�nite element approaches: Part 1. Use of cubic splines in �nding minimal surface membranes. International journal for numerical methods in engineering
56, 5 (2003), 651–668.
Oleksiy Busaryev, Tamal K. Dey, Huamin Wang, and Zhong Ren. 2012. Animating
bubble interactions in a liquid foam. ACM Trans. Graph. 31, 4 (2012), 63.
Fang Da, Christopher Batty, and Eitan Grinspun. 2014. Multimaterial Mesh-Based
Surface Tracking. ACM Trans. on Graphics (SIGGRAPH 2014) 33 (2014), 1–11.
Fang Da, Christopher Batty, Chris Wojtan, and Eitan Grinspun. 2015. Double Bubbles
Sans Toil and Trouble: Discrete Circulation-Preserving Vortex Sheets for Soap Films
and Foams. ACM Trans. on Graphics (SIGGRAPH 2015) 34, 4 (2015), 149:1–149:9.
Fang Da, Christopher Batty, Chris Wojtan, and Eitan Grinspun. 2017. Project pageof "Double Bubbles Sans Toil and Trouble: Discrete Circulation-Preserving Vortex Sheets for Soap Films and Foams". (2017). http://www.cs.columbia.edu/cg/
doublebubbles/
Fang Da, Christopher Batty, Chris Wotjan, and Eitan Grinspun. 2015. Use of Fast
Multipole to Accelerate Discrete Circulation-Preserving Vortex Sheets for Soap
Films and Foams. (2015), 2.
Mathieu Desbrun, Mark Meyer, Peter Schröder, and Alan H. Barr. 1999. Implicit Fairingof Irregular Meshes Using Di�usion and Curvature Flow. In Proceedings of the 26th
Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH ’99).
ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, 317–324.
Kong Dexing, Liu Kefeng, and Wang Zenggui. 2009. Hyperbolic mean curvature �ow:
evolution of plane curves. Acta Mathematica Scientia 29, 3 (2009), 493–514.
Jesse Douglas. 1931. Solution of the problem of Plateau. Trans. Amer. Math. Soc. 33, 1(1931), 263–321.
Roman Durikovic. 2001. Animation of Soap Bubble Dynamics, Cluster Formation and
Collision. Comput. Graph. Forum 20, 3 (2001), 67–76.
Raymond Yun Fei. 2017. Personal communication. (Aug. 2017).
Luca Giomi and Lakshminarayanan Mahadevan. 2012. Minimal surfaces boundedby elastic lines. Proceedings of the Royal Society A: Mathematical, Physical andEngineering Sciences 468, 2143 (mar 2012), 1851–1864.
Giulio G Giusteri, Luca Lussardi, and Eliot Fried. 2017. Solution of the Kirchho�–Plateau
problem. Journal of Nonlinear Science 27, 3 (2017), 1043–1063.
Gaël Guennebaud, Benoît Jacob, et al. 2010. Eigen v3. http://eigen.tuxfamily.org. (2010).
Jenny Harrison. 2014. Soap �lm solutions to Plateau’s problem. Journal of GeometricAnalysis 24, 1 (2014), 271–297.
Jenny. Harrison and Harrison Pugh. 2015. Plateau’s Problem: What’s Next. ArXiv
e-prints (Sept. 2015), 27. arXiv:1509.03797
Chun-Lei He, De-Xing Kong, and Kefeng Liu. 2009. Hyperbolic mean curvature �ow.
Journal of Di�erential Equations 246, 1 (2009), 373–390.
Jeong-Mo Hong, Ho-Young Lee, Jong-Chul Yoon, and Chang-Hun Kim. 2008. Bubbles
alive. ACM Trans. Graph. 27, 3 (2008), 48:1–48:4.
Michael Hutchings, Frank Morgan, Manuel Ritoré, and Antonio Ros. 2002. Proof of the
double bubble conjecture. Annals of Mathematics 155, 2 (2002), 459–489.
Cyril Isenberg. 1978. The science of soap �lms and soap bubbles. Courier Corporation,
North Chelmsford, Massachusetts, USA.
Alec Jacobson, Daniele Panozzo, et al. 2013. libigl: A simple C++ geometry processing
library. (2013). http://igl.ethz.ch/projects/libigl/.
Wenzel Jakob. 2010. Mitsuba renderer. (2010). http://www.mitsuba-renderer.org.
Byungmoon Kim, Yingjie Liu, Ignacio Llamas, Xiangmin Jiao, and Jarek Rossignac. 2007.
Simulation of Bubbles in Foam With The Volume Control Method. ACM Trans. on
Graphics (SIGGRAPH 2007) 26, 3 (2007), 10.
Erwin Kreyszig. 1991. Di�erential geometry. Dover Publications, New York.
Philippe G LeFloch and Knut Smoczyk. 2008. The hyperbolic mean curvature �ow.
Journal de mathématiques pures et appliquées 90, 6 (2008), 591–614.
Rafael López. 2013. Constant mean curvature surfaces with boundary. Springer Science& Business Media, New York, NY, USA.
Carlo Mantegazza. 2011. Lecture notes on mean curvature �ow. Vol. 290. Springer
Science & Business Media, New York, NY, USA.
Barry Merriman, James K. Bence, and Stanley J. Osher. 1994. Motion of Multiple
Junctions: A Level Set Approach. J. Comput. Phys. 112 (1994), 334–363.
Mark Meyer, Mathieu Desbrun, Peter Schröder, and Alan H. Barr. 2002. DiscreteDi�erential-Geometry Operators for Triangulated 2-Manifolds. (2002).
Marek Krzysztof Misztal, Kenny Erleben, Adam Bargteil, Jens Fursund, Brian Bunch
Christensen, Jakob Andreas Bærentzen, and Robert Bridson. 2014. Multiphase�ow of immiscible �uids on unstructured moving meshes. IEEE transactions on
visualization and computer graphics 20, 1 (2014), 4–16.
Rhudaina Z. Mohammad and Karel Švadlenka. 2002. Proof of the double bubble conjecture. Ann. Math (2002), 0406017.
Matthias Müller. 2009. Fast and robust tracking of �uid surfaces. In Proceedings of the2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA ’09).
ACM, New York, NY, USA, 237–245.

199:11

Hao Pan, Yi-King Choi, Yang Liu, Wenchao Hu, Qiang Du, Konrad Polthier, Caiming
Zhang, and Wenping Wang. 2012. Robust Modeling of Constant Mean Curvature
Surfaces. ACM Trans. Graph. 31, 4, Article 85 (July 2012), 11 pages.
Saket Patkar, Mridul Aanjaneya, Dmitriy Karpman, and Ronald Fedkiw. 2013. A hybrid
Lagrangian-Eulerian formulation for bubble generation and dynamics. In The ACM
SIGGRAPH / Eurographics Symposium on Computer Animation, SCA ’13, Anaheim,
CA, USA, July 19-21, 2013. ACM Press, New York, NY, USA, 105–114.
Jesus Perez, Miguel A. Otaduy, and Bernhard Thomaszewski. 2017. Computational
Design and Automated Fabrication of Kirchho�-Plateau Surfaces. ACM Trans. on
Graphics (Proc. of ACM SIGGRAPH) 36, 4 (2017), 62.1–62.12.
Ulrich Pinkall and Konrad Polthier. 1993. Computing Discrete Minimal Surfaces andTheir Conjugates. Experimental Mathematics 2 (1993), 15–36.
Konrad Polthier and Wayne Rossman. 2002. Discrete Constant Mean Curvature Surfaces
And Their Index. In Visualization and Mathematics. Springer Verlag, 47–77.
Thomas R Powers, Greg Huber, and Raymond E Goldstein. 2002. Fluid-membrane
tethers: minimal surfaces and elastic boundary layers. Physical Review E 65, 4 (2002),
041901.
Giuseppe Pucci, Daniel M Harris, and John WM Bush. 2015. Partial coalescence of soap
bubbles. Physics of Fluids 27, 6 (2015), 061704.
Tibor Rado. 1930. On Plateau's Problem. The Annals of Mathematics 31, 3 (jul 1930),
457.
Robert I. Saye and James A. Sethian. 2013. Multiscale Modeling of Membrane Rearrangement, Drainage, and Rupture in Evolving Foams. Science 340, 6133 (2013),
720–724.
Jos Stam. 1999. Stable Fluids. In Proceedings of the 26th Annual Conference on Computer
Graphics and Interactive Techniques (SIGGRAPH ’99). ACM Press/Addison-Wesley
Publishing Co., New York, NY, USA, 121–128.
Michael Struwe. 2014. Plateau’s Problem and the Calculus of Variations.(MN-35). Princeton University Press, Princeton, New Jersey, USA.
John M Sullivan and Frank Morgan. 1996. Open problems in soap bubble geometry.
International Journal of Mathematics 7, 06 (1996), 833–842.
Jean E Taylor. 1976. The structure of singularities in soap-bubble-like and soap-�lm-likeminimal surfaces. Annals of Mathematics 103, 3 (1976), 489–539.
Nils Thürey, Chris Wojtan, Markus Gross, and Greg Turk. 2010. A multiscale approach
to mesh-based surface tension �ows. ACM Transactions on Graphics 29, 4 (2010), 10.
Karel Švadlenka, Elliott Ginder, and Seiro Omata. 2014. A variational method for
multiphase volume-preserving interface motions. J. Comput. Appl. Math. 257 (2014),
157–179.
Shing-Tung Yau. 2000. Review of geometry and analysis. Mathematics: frontiers andperspectives (2000), 353–401.
Yizhong Zhang, Huamin Wang, Shuai Wang, Yiying Tong, and Kun Zhou. 2012. A
deformable surface model for real-time water drop animation. IEEE Transactions on
Visualization and Computer Graphics 18, 8 (2012), 1281–1289.
Wen Zheng, Jun-Hai Yong, and Jean-Claude Paul. 2006. Simulation of Bubbles. In
Proceedings of the 2006 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (SCA ’06). Eurographics Association, Aire-la-Ville, Switzerland, Switzerland,
325–333.
Bo Zhu, Ed Quigley, Matthew Cong, Justin Solomon, and Ronald Fedkiw. 2014. Codimensional surface tension �ow on simplicial complexes. ACM Transactions on
Graphics (TOG) 33, 4 (2014), 111.

ACM Transactions on Graphics, Vol. 36, No. 6, Article 199. Publication date: November 2017.