Continuum Foam:
A Material Point Method for Shear-Dependent FlowsYONGHAO YUE and BREANNAN SMITH
Columbia University

CHRISTOPHER BATTY
University of Waterloo

CHANGXI ZHENG and EITAN GRINSPUN
Columbia University
We consider the simulation of dense foams composed of microscopic bubbles, such as shaving cream and whipped cream. We represent foam not
as a collection of discrete bubbles, but instead as a continuum. We employ
the Material Point Method (MPM) to discretize a hyperelastic constitutive
relation augmented with the Herschel-Bulkley model of non-Newtonianviscoplastic flow, which is known to closely approximate foam behavior.
Since large shearing flows in foam can produce poor distributions of material
points, a typical MPM implementation can produce non-physical internalholes in the continuum. To address these artifacts, we introduce a particle
resampling method for MPM. In addition, we introduce an explicit tearing
model to prevent regions from shearing into artificially-thin, honey-likethreads. We evaluate our method’s efficacy by simulating a number of dense
foams, and we validate our method by comparing to real-world footage of
foam.
Categories and Subject Descriptors: I.3.7 [Computer Graphics]: ThreeDimensional Graphics and Realism—Animation; I.6.8 [Simulation andModeling]: Types of Simulation—Animation
General Terms: Algorithms
Additional Key Words and Phrases: Foam, material point method, particle
resampling, tearing, shear thinning, shear thickening, viscoplasticity

1.

INTRODUCTION

Foams are a unique class of materials that arise in many problem domains and that exhibit a remarkable range of behaviors. Application
domains as diverse as firefighting, oil recovery, chemical filtration,and even industrial textile processing employ foams [Weaire and
Hutzler 2001]. In computer graphics, foams are a useful componentof many animated scenarios; indeed, the recent Pixar animated short“Partysaurus Rex” prominently showcases soap bubbles and foams[Harder and Mangnall 2013]. The simulation of foam is thus an
important problem in both engineering and graphics contexts.
Foam bubbles encompass a wide spectrum of spatial scales, ranging from the macroscopic bubbles formed by bath soaps on the orderof centimeters, to the immense number of bubbles that compose
shaving foams on the order of tens of micrometers. We focus on therelatively unexplored case of dense foams, exemplified by materials
such as shaving cream and whipped cream. As our focus is on dense
foams, we do not model the drainage [Koehler et al. 1998] andcollapse [Garrett 1993] effects characteristic of macro-scale bubbles.
To guide the design of our method, we first consider the empiricalproperties of foams.
Dense foams consist of myriad individual bubbles, each individually difficult to see due to its small size; a typical shaving cream
bubble is 40µm in diameter [Ovarlez et al. 2010]. Despite consist-ing of many small bubbles interacting with each other, the observed
macroscopic behavior of a dense foam is that of a seamless continuum. This observation suggests that we can simulate foam in a
homogenized way to resolve large-scale and visually compelling
dynamics.
When viewed as a continuous medium, foam exhibits a nonNewtonian force response. Due to the microscopic stochastic rearrangement of bubbles, foams exhibit viscoplasticity, meaning they
undergo rate-dependent permanent deformation. Foams also tend to
exhibit shear thinning effects, in which the continuum flows moreeasily under larger applied stresses. Shear thinning, which also occurs in blood, motor oil, and paint, plays a significant role in theutility of many materials. For example, shear thinning allows paintto readily flow off of a brush, but not off of a wall. Viscoplasticity
and shear thinning can induce large shearing flows and topological changes with the potential to radically rearrange a material’s
configuration.
With these observations in mind, we simulate bulk foam as acontinuous medium, allowing us to treat larger volumes of foamthan would be tractable at the individual bubble level. We adopt the
Herschel-Bulkley constitutive relation for foam [Weaire and Hutzler 2001]. While a precise understanding of the relation between
microscopic bubble mechanics and the aggregate behavior of foamat the macroscopic scale remains an open problem [Weaire 2008],
empirical experiments have shown the Herschel-Bulkley model tobe a highly effective phenomenological approximation. As an added
bonus, we will demonstrate that the Herschel-Bulkley constitutive relation is general enough to model a wide variety of other interestingmaterials, including shear thickening “oobleck.”
To address the challenges posed by foam’s extreme shape changes,
we adopt the hybrid particle/grid method known as the Material
Point Method (MPM), which has been noted for the relative easewith which it handles large deformations and topology changes[Stomakhin et al. 2013; Stomakhin et al. 2014]. A direct applicationof MPM to our setting falls short, however. Large shearing flows
can produce poor material point distributions, creating non-physicalinternal holes (see Figure 6). We remedy this problem with a novel
particle resampling method for MPM that maintains uniform particle distributions while retaining the physical properties and thegeometry of the foam.
Finally, we observe that an explicit tearing model is required
to realistically simulate foam; without a tearing model, regions offoam can become progressively thinner, yielding artificially-slender,
honey-like threads. These unnatural threads form because the grid
maintains artificial topological connections at the grid resolution
through the stencil of the shape functions; the particles are only
regarded as topologically distinct when they are separated by a
ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

2•

Y. Yue et al.
thaler and Pajarola 2008], though they have not been applied tofoam structures.
Bubble Films as Triangle Meshes. An alternate approach is tot = 0.63s

t = 0.63s

t = 0.90s

t = 0.90s

model bubble interfaces as dynamic triangle meshes. This is exemplified by the Surface Evolver model [Brakke 1992], which seeks
equilibrium configurations of surfaces under surface tension andvarious constraints including volume preservation. Ďurikovič [2001]
developed a different triangle surface model for dynamic bubbles,
treating the film as an elastic deformable object with surface tension
and using Lennard-Jones forces for bubble interactions.
Particle Models for Foams and Bubbles. To scale up to denser

t = 1.08s

t = 1.08s

t = 1.50s

t = 1.50s

Fig. 1. Comparison of our method to real-world video footage of adollop of shaving cream falling from an oscillating block. Our simulation
results (left of each pair with blue labels) match closely with the captured
video footage (right of each pair with pink labels).

sufficiently large distance. We therefore propose an extension to
directly detect and tear weak regions based on their accumulated
plasticity.
In this work we offer three core contributions: we extend the MPM
framework to support the general viscoplastic Herschel-Bulkley material model, we introduce an effective particle resampling method
for MPM, and we introduce an explicit tearing model. Taken together
these features enable our continuum-based method to realisticallysimulate the familiar behaviors of dense foams, as well as a wide
variety of related materials. We demonstrate the method’s efficacy
with several animated examples and with a comparison to recordingsof real foam (see Figure 1).

2.

RELATED WORKThe Physics of Foam. Due to the wide array of domains in whichfoams play a role, there has been substantial research into their properties and behavior. The text by Weaire and Hutzler [2001] provides
a broad overview of earlier theoretical, experimental, and computational approaches to foam mechanics. The field of rheology is
the study of flowing matter, including foams, although a systematicrheological understanding of foams remains to be fully developed[Weaire 2008]. As Weaire notes, the dominant continuum model isthe nonlinear Herschel-Bulkley fluid model, a generalized modelfor non-Newtonian fluids that we adopt. While the aggregate behavior of foam is ultimately dictated by the small-scale mechanics of
its constituent bubbles, for many applications, including ours, theHerschel-Bulkley model is well-suited.
Navier-Stokes Approaches. Multiphase fluid simulation meth-ods [Hong and Kim 2005; Losasso et al. 2006] have been adapted tosimulate the detailed deformations of small collections of bubbles[Zheng et al. 2006; Kim et al. 2007; Kim et al. 2010]. Multiphaseextensions of the level set method are used to distinguish individual
pockets of air, with the Navier-Stokes equations applied to evolve thecomplete system. Sethian et al. [2013] recently extended a similar
framework with more advanced drainage and film rupture modeling.
Multiphase Smoothed Particle Hydrodynamics (SPH) techniques
can likewise model submerged bubbles [Müller et al. 2005; Solen-

foams composed of small bubbles, for which the deformations ofindividual films become less relevant, many authors have adopted
models in which each bubble is represented by a single particle,
typically with spring-like interparticle forces. These approaches areoften layered atop a single-phase fluid solver in order to increase theapparent level of detail and realism. Focusing on foams themselves,Kück et al. [2002] applied spring forces between bubble particles
to approximate the Plateau conditions of soap films at equilibrium.
Greenwood and House [2004] later coupled this model to a NavierStokes fluid solver, while Thürey et al. [2007] coupled particle-based
bubbles to a shallow water model. Cleary et al. [2007] augmented an
SPH liquid simulator with a discrete bubble model, using cohesion
forces to encourage the formation of foam “rafts” on the liquidsurface, while Ihmsen et al. proposed a method to layer foam onto
an SPH animation as a post-process [Ihmsen et al. 2012]. To better
approximate the foam geometry and to improve volume preservation,
Busaryev et al. [2012] exploited the properties of Voronoi diagrams.
Hong et al. coupled an SPH bubble model with an Eulerian liquid toenliven the bubble behavior [Hong et al. 2008]. Patkar et al. [2013]
tightly coupled a particle-based subgrid compressible model to a
two-phase fluid simulator.
Multiphase Continua and Homogenization. Similar in spirit to

our work, Kim et al. [2010] took a step away from individual interparticle forces in the context of dispersed bubble flows, instead
adopting a variable density multiphase continuum model that accounts for combined bubble/liquid behavior. Nielsen et al. [2013]
applied a continuum model to a similar problem in which small
droplets of water are dispersed in air to produce a spray.
Homogenization methods [Nesme et al. 2009; Kharevych et al.
2009] develop coarse-scale discretizations that approximate aggregate fine-scale inhomogeneities in material properties and geometry.
The extreme nonlinearities exhibited by foams makes this strategy
inapplicable.
Methods for Viscoplastic and Non-Newtonian Materials. Var-ious non-Newtonian and viscoplastic models have been developed
to animate materials whose behavior cannot be neatly categorized
as fluid or solid. Early on, Terzopoulos and Fleischer [1988] developed deformable models that incorporated viscoplasticity andfracturing effects. Goktekin et al. [2004] gradually accumulated
strain within an Eulerian fluid solver to add elastic forces, while
Bargteil et al. [2007] augmented a Lagrangian finite element method(FEM) elasticity solver with remeshing to robustly capture dramatic
viscoplastic flows. Subsequent work in the Lagrangian FEM setting
has focused on accelerating remeshing and reducing the associatednumerical smearing [Wojtan and Turk 2008; Wicke et al. 2010].
Point-based methods using both SPH and moving-least-squares(MLS) approaches have also been extended to non-Newtonian materials [Müller et al. 2004; Keiser et al. 2005; Solenthaler et al. 2007;
Gerszewski et al. 2009].

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Continuum Foam: A Material Point Method for Shear-Dependent Flows•3

Independent of the particular underlying numerical technique,graphics authors have adapted a range of constitutive models for capturing non-Newtonian and viscoplastic effects. O’Brien et al. [2002]
proposed a simple additive plasticity model, while later authors
adopted a more powerful multiplicative plasticity model [Irving
et al. 2004]. Bargteil et al. [2007] added support for creep effects
and work hardening and softening. While we are motivated by foamin particular, the constitutive model we propose offers increased
flexibility to handle a variety of plastic effects, including perfectplasticity, viscoplasticity, shear thinning and shear thickening.

“permanent” fracture, and can be overkill for foam simulation. To
this end, we propose a simple technique to explicitly model tearing
in foams that also allows foams to reconnect.

FLIP and MPM Methods. We build on a hybrid particle/gridmain Ω ⊂ R3 , with scalar density field ρ(x, y, z), vectorial velocity
field v(x, y, z), tensorial Cauchy stress field σ(x, y, z), and subject
to gravitational acceleration g. The dynamics of this medium evolve
over time t under the Euler-Lagrange partial differential equation

framework for elastoplastic materials called the Material Point
Method (MPM) [Sulsky et al. 1995]. MPM uses Lagrangian particles
to handle the transport of information through space, but exploits
a background regular grid structure to compute elastic forces. As
we discuss below (§3), this confers many of the advantages of both
Lagrangian and Eulerian methods. MPM was originally developed
as an application of the related FLIP method for compressible flow
to elastic media [Brackbill and Ruppel 1986]. As with FEM, MPM
is derived in the weak form and can treat history-dependent quantities as material points. Zhu and Bridson applied FLIP to reduce
dissipation in incompressible fluid animations [Zhu and Bridson2005]; notably, they also abstracted away individual grains to apply
a continuum model for granular media, later extended by Narain etal. [2010]. Recently, Stomakhin et al. [2013] developed an MPM
technique and constitutive model to animate the unique behavior ofsnow, likewise treating it as a continuum. Stomakhin et al. [2014]
subsequently extended their MPM technique to incorporate phasechanges due to heat flow and to enforce incompressibility with aprojection.
Sampling. Nonuniform point distributions can pose problems fornumerical methods that advect particles. Closest to our work is thatof Edwards et al. [2012], who pointed out that shearing flows canquickly introduce anisotropic point distributions for particle-in-cell
methods, like MPM; they proposed a (non-conservative) method toresample points. Stomakhin et al. [2014] also pointed to the needfor a resampling technique for MPM as a direction for future work.
Particle-based vortex methods [Cottet and Koumoutsakos 2000] canalso suffer from this problem. Particle position adjustment, particle
resampling, and the use of source/sink terms have all been suggested
in the context of FLIP methods [Losasso et al. 2008; Narain et al.
2010; Ando et al. 2012]. Surface representation methods, such as
particle/level set hybrids [Enright et al. 2002; Hieber and Koumoutsakos 2005], and surfel-based methods [Keiser et al. 2005], canalso require resampling to avoid holes or loss of accuracy. A largely
orthogonal but related problem arises in adaptive simulation, forboth FLIP and SPH, in which particles must be merged or split toachieve different scales (e.g., [Adams et al. 2007; Ando et al. 2012]).
Schechter and Bridson [2012] proposed using Poisson-disk sampling to place ghost air particles at each step and to seed the initial
SPH particles according to a level set description of the geometry;
we build on these ideas to resample the interior region of our MPM
foam.
Material Failure and Fracture. Foams exhibit behaviors charac-

teristic of both fluids and solids, leading to phenomena such as theformation of stiff peaks. Similar to solids, foams can tear and form
stiff peaks, yet similar to fluids, detached samples of foam can reconnect under contact. While methods to explicitly handle material
failure and fracture for stiff materials have been proposed for MPM[Zhou 1998; Schreyer et al. 2002; Sulsky and Schreyer 2004; Banerjee 2004; Guo and Nairn 2006], these techniques typically consider

3.

BACKGROUND: MATERIAL POINT METHOD

Notation. In the following, we depict scalars in regular romantype (e.g., ρ and J), vectors in bold italic type (e.g., v and g), andrank-two tensors in bold sans serif type (e.g., σ and s).
Continuous Setting. Consider a continuous medium over a do-Dv= ∇ · σ + ρg
Dt
subject to the mass conservation conditionρDρ+ ρ∇ · v = 0
DtD·
where Dt
denotes the material derivative, i.e., Dφ= ∂φ+ v · ∇φDt∂t
for a time-varying field φ(x, y, z). The details of our constitutive
model, which defines σ, are presented in §5.

Discretization. The material point method [Sulsky et al. 1995]

discretizes the Euler-Lagrange equations in both space and time
via a sequence of Eulerian and Lagrangian steps. This discretization approach is ideally suited for non-Newtonian and viscoplastic
materials, which possess characteristics of both fluids and solids.
Material geometry and kinematic properties (density, velocity, deformation gradient, etc.) are associated with particles, which allows
them to be transported easily through space in a Lagrangian manner. This circumvents traditional challenges of Eulerian methods,
especially artificial numerical dissipation.
Material forces and particle accelerations are computed using astatic regular (Cartesian) background grid, and then applied to the
Lagrangian particles. Recourse to a static grid avoids the potentiallycomplex and expensive remeshing required by methods that transport the mesh/grid itself in a Lagrangian manner, while the regulargrid structure simplifies discretization and optimization. The useof a grid to increment rather than overwrite particle data reduces
numerical dissipation.
MPM also makes various geometric aspects of large-deformation
foam dynamics easier to handle: the Lagrangian particles implicitly
account for topological splitting and merging; the Eulerian grid
implicitly accounts for (self-)collision.
The net result is a method that is highly effective at modelingboth solid and fluid properties. Like mesh-based Lagrangian techniques, MPM can be used with an arbitrary constitutive law tomodel complex materials; this flexibility allowed Stomakhin et al.
[2013] to achieve impressive snow effects by combining MPM witha specially-designed snow constitutive model. Like many discretization methods, MPM is also able to provide exact conservation ofcertain quantities. In our implementation, which generally follows
Stomakhin et al. [2013], mass is exactly conserved.

4.

OVERVIEW

We begin by summarizing the procedures that comprise each time
step of our method. Figure 2 depicts this sequence, which mirrors
the overall structure of the method of Stomakhin et al. [2013]. We
make substantial modifications in the procedures highlighted in

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

4•

Y. Yue et al.
Particle domain

Particle state
2. Tear detection (§7.1)

7. Deformation gradient (§5.2)

8. Deformation gradient limiting (§7.2)

9. Particle velocities

1. Grid velocities and mass

(11) Compute plastic recovery. We relax the plasticity history to
account for the strengthening of bonds between recently adhered bubbles (§7.3).

3. Grid forces (§5)(12) Remove subgrid geometry. We remove thin, weakened particles whose physics cannot be resolved on the grid (§7.4).

Grid domain

4. Force limiting (§7.2)

5. Grid velocities
6. Grid collisions

10. Plastic flow (§5.2)

11. Plastic recovery (§7.3)

12. Subgrid geometry removal (§7.4)
13. Particle collisions

14. Position update

15. Avoid a void (§6)

Fig. 2. An overview of the stages of the MPM algorithm, with our
additions highlighted in pink. The arrows indicate the data flow.

pink, which implement our constitutive model (procedures 3 and10), our tearing model (procedures 2, 4, 8, 11, 12), and our particle
resampling technique (procedure 15).
Each particle carries several pieces of data: position x, velocityv, mass m, (initial) volume V , deformation gradient F, and noremalized elastic Green strain tensor b = det(be )−1/3 be . Here andhenceforth the bar operator isotropically rescales a tensor to haveunit determinant, i.e., x = det(x)−1/3 x. Our Eulerian grid-based
FEM discretization relies on the one-dimensional cubic B-spline
basis functions advocated by Stomakhin et al. [2013].
Briefly, the steps of our algorithm are:(1) Rasterize particle data to the grid. At the start of every time
step we transfer each particle’s mass and velocity to the Eulerian
grid.
(2) Detect tearing regions. We examine each particle’s accumulated plasticity to detect regions that are tearing; this determines
where the material should be weakened during the next step(§7.1).
(3) Compute grid forces. Returning to the grid, we compute thegrid-based forces according to our constitutive model (§5).
(4) Limit forces. We modify the stress and the force in tearing
regions according to our tearing model (§7.2).
(5) Update grid velocities. We integrate the resulting grid-based
forces to yield new grid-based velocities.
(6) Resolve grid-based collisions. We update the grid-based velocities to account for collisions with objects ([Stomakhin et al.
2013] §8).
(7) Update particle deformation data. We update the per-particle
deformations based on the velocities (§5.2).

(8) Limit deformation gradient. We limit the change in the deformation gradient in the tearing regions according to our tearing
model (§7.2).
(9) Update particle velocities. We apply the incremental change
in the grid velocities back to the particle velocities (à la FLIP
methods).
(10) Compute plastic flow. We compute a plastic update fromthe deformation gradient, transferring excess elastic strain into
permanent plastic strain (§5.2).

(13) Resolve particle-based collisions. We update the particlebased velocities to account for collisions with objects ([Stomakhin et al. 2013] §8).
(14) Update particle positions. We update the particle positions
according to the particle velocities.

(15) Avoid a void. Finally, we resample the material points to fill
non-physical voids and to maintain a uniform particle distribution (§6).

5.

CONSTITUTIVE MODEL

Foams (such as Figure 3, bottom) consist of concentrated dispersionsof gas bubbles immersed in a surfactant solution [Weaire and Hutzler
2001]. Foams exhibit features over a broad range of length scales,from the microscopic liquid-gas interfaces of individual bubbles, tothe macroscopic body which is well approximated as a continuous
medium.
At the macroscopic scale, foams behave as shear thinning fluids
for which the relationship between shear stress and strain rate is wellapproximated by the Herschel-Bulkley phenomenological law [Herschel and Bulkley 1926; Cohen-Addad et al. 2013]. As described by
Weaire and Hutzler [2001], the application of a small stress deforms
individual bubbles, which is reflected at a macroscopic scale as areversible, elastic behavior. Increasing the stress, however, induces
a rearrangement of bubbles, which is in turn reflected as a macroscopic plastic deformation, or a change in rest state. A continued
increase in stress eventually induces a shear thinning flow of thesurfactant solution and its immersed bubbles.
We model foam macroscopically. To capture a wide range of
observed foam phenomena, we desire a constitutive relation that canmodel elasticity, plasticity and shear thinning flow.
Deformation Gradient and Tensor. We begin our presentationof the constitutive model with the deformation gradient, F(x), a
linear mapping of an infinitesimal material neighborhood from its
undeformed to its deformed position [Simo and Hughes 1998]. The
deformation gradient maps every infinitesimal vector δx anchoredat material point x in the undeformed configuration to its corresponding vector F(x)δx in the deformed configuration. We typically omit the argument, x = (x, y, z), while bearing in mind thatF : Ω 7→ R3×3 is in general a heterogeneous field over the material
domain.
In the treatment of plasticity, it is convenient to decompose F into
parts associated to the elastic, Fe , and the plastic, Fp , deformation,
via the (multiplicative) decomposition [Simo and Hughes 1998;
Irving et al. 2004; Bargteil et al. 2007; Wicke et al. 2010; Jones et al.
2014]
F = Fe Fp .
Plasticity serves as the macroscopic description of a volumepreserving bubble rearrangement, thus by ansatz det(Fp ) = 1.
In addition, our treatment of elasticity will also refer to the left
Cauchy-Green tensor b = FFT and its elastic part be = Fe Fe T .
The tensors b and be are objective (i.e., frame-indifferent).
Elasticity. The application of a slight strain on the foam performs

reversible elastic work, captured by the hyperelastic stored energy

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

(1)whose terms account for volumetric and shear deformation,weighted by the material-dependent bulk modulus
κ and the shear
R
modulus µ, respectively; correspondingly, Ω W dVol is the stored
energy, which is also objective.
The volume-dependent energy density11
Wv (J) = κ (J 2 − 1) − ln J
22
is expressed in terms of the square and natural logarithm of thescalar J = det(F), whereas the shear-dependent energy density1eµ(Tr[b ] − 3)2is expressed in terms of the volume preserving left-Cauchy Greenetensor b = det(be )−1/3 be = J −2/3 be . This energy function
alone characterizes the constitutive relation for finite stresses. In
particular, it gives rise to the Kirchhoff stress tensor (see Appendix B
for the derivation)eWs (b ) =

∂W eT
κ
e
(2)
F = (J 2 − 1)I + µ dev[b ] ,∂Fe2whose two terms encode an isotropic response resulting from volumetric and shear strains, respectively. In the expression above,dev[x] = x − Tr[x]I, and I denotes the 3 × 3 identity tensor.
3
For any tensor x, the deviatoric operator dev[x] returns the tracefree part of x, i.e., Tr[dev[x]] = 0. Thus, dev[I] = 0, anddev[dev[x]] = dev[x]. We will use all of these identities in §5.2.
The Kirchhoff tensor is related to the Cauchy stress by σ = τ /J.
The internal force f is given by f = ∇ · σ.
The treatment of plasticity shall make use only of the deviatoricepart of the Kirchhoff stress tensor, s = dev[τ ] = µ dev[b ], its(scalar) magnitude s = kskF , and its normalized form ŝ = s/s;
here k · kF denotes the Frobenius norm. Intuitively, s encodes themagnitude of the material’s elastic response to shear, ignoring theorientation of the shear mode.
τ =Plastic Yield. The von Mises yield condition [Simo and Hughes
1998; O’Brien et al. 2002] makes precise the limits of the elastic
regime, in terms of a material-dependent yield stress σY ; note that
because we employ an isotropic elastic response, σY is a scalar. So
long as the conditionr
2
σY ≤ 0(3)
Φ(s) = s −3

is satisfied, the material response remains elastic. Intuitively, thematerial cannot sustain a purely elastic response when the magnitudeof the elastic response to shear deformation exceeds the thresholdσY . We refer to the quantity max(0, Φ(s)) as the yield excess, andthe inequality Φ(s) ≤ 0 as the yield condition.
For more general plastic behavior, work hardening and softening
are modeled by an update of the yield stress [Simo and Hughes
1998]; however, in quasistatic experiments of foams [Weaire and
Hutzler 2001] it was observed that there is no apparent tendency for
hardening or softening. We therefore safely neglect this effect.

5.1

Plastic Model

When the von Mises condition is violated, the foam deforms plastically; mathematically, this is modeled by updating the plastic strain

simulation resultse

W = Wv (J) + Ws (b )5Real footage Herschel-Bulkley viscoplastic

density ([Simo and Hughes 1998]-§9.2)•

perfect plastic

Continuum Foam: A Material Point Method for Shear-Dependent Flowst = 0.50s

t = 0.63s

t = 0.77s

t = 0.90s

t = 1.08s

Fig. 3. Comparison of real footage (bottom row) of a dollop of dense
shaving cream falling from an oscillating block to simulations with threeplasticity models (perfect plastic, viscoplastic and Herschel-Bulkley in the
top three rows). The behavior of the Herschel-Bulkley model most closely
matches that of the real experiment. For the sake of comparison, we only
vary the yield stress σY , the viscosity η, and the Herschel-Bulkley power h.
The thin filament visible in the perfect plasticity model can be removed bymodifying the tearing parameters.

(and thus reducing the elastic strain) in the decomposition of thestrain tensor. The time derivative of be is given by the identityḃe = Lbe + be LT + Lv be ,(4)

2
Lv be = − Tr[be ]γ ŝ ,3(5)which involves both the spatial velocity gradient L = ∇v =ḞF−1 , as well as the Lie derivative of be , Lv be (see Simo andHughes [1998]). The former captures the change in the strain due
to the flow field itself, while the latter Lie derivative captures thechange relative to the flow induced by v, i.e., due to the plastic flow.
In rheology, one adopts a flow rule to describe the plastic deformation, in which excess elastic strain flows into plastic strain over time.
We use the flow rule of Simo and Hughes (see [Simo 1988]-§1.4,[Simo and Hughes 1998]-§9)

where γ is the flow rate, and ŝ is the flow direction. This flow rule isobjective and can be derived by applying the principle of maximum
plastic dissipation to the stored energy function (1) and the yield
condition (3) [Simo 1988]-§1.4.
Perhaps the simplest update model is perfect plasticity, in which
the excess elastic strain (i.e., the amount that exceeds the von Mises
condition) is instantaneously transferred from the elastic to theplastic strain; in the context of perfect plasticity, γ in Eq. (5) iscalled the consistency parameter, and is obtained by solving Eq. (5)with the condition that the elastic strain, and the attendant stress,
never exceeds the yield stress (see Figure 4).
Albeit simple, perfect plasticity fails to capture the viscoplastic
behavior of foam. Consider an experiment in which we place real

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

6•

Y. Yue et al.

stressstressyield stressyield stress

Viscoplastic

effective viscosity[Bargteil et al. 2007]
Viscoplastic

Perfect plastic

Perfect plastic

HerschelBulkley
stressstraintime

Fig. 4. Comparison of constitutive relations: (Left) An illustration of thestress-strain relation for the perfectly plastic model. (Right) A comparison of
the perfect plastic and viscoplastic models as we increase the total strain at aconstant strain rate. While the perfectly plastic model immediately removes
any excess stress, the viscoplastic model acts over a finite timescale. Hence,
under the viscoplastic model, the total stress can exceed the yield stress.
This (large) stress is subject to elastic deformation, hence the material resists
flowing.

foam on the underside of a wooden block (see Figure 3 and theaccompanying video) and shake vigorously in a horizontal motion.
Initially, the foam vibrates in an elastic manner, and the material
flows slowly; over time, however, the foam begins to rapidly flow,
leading to the formation of a thin neck that eventually tears away.
Perfect plasticity fails to capture this behavior: the foam flows fromthe onset of agitation.
The block experiment motivates the need for a viscoplastic model(see Figure 4), where the flow rate is replaced by a function of theviscosity-normalized yield excess asΦ(s)γ = γ(s) = max 0,,(6)ηwhere η is the viscosity (or relaxation) coefficient. Eq. (6) is known
as the Bingham plastic model [Bingham 1922], a common modelfor toothpaste, mayonnaise, and similar materials.
With a viscoplastic model in hand, we revisit the block experiment(see Figure 3). When the foam is agitated, the excess strain cannot be
instantaneously eliminated, and thus early vibrations are seemingly
elastic; over a time scale determined by the viscosity η, the ratelimited plastic flow leads to the formation of a neck and then totearing.
While a considerable improvement over perfect plasticity, a viscoplastic law does not capture the shear thinning behavior of foam;
the formation of the neck happens too slowly. Refer to Fig. 5 andthe accompanying video for a comparison of the plasticity models.
Herschel-Bulkley Model. The Herschel-Bulkley model re-places (6) with a power law relating the plastic flow rate to theyield excess:
 
1/h
Φ(s)
γ(s) = max 0,η(7)Shear thinning and thickening materials are characterized by h < 1
and h > 1, respectively.
To reveal the advantage of the Herschel-Bulkley model, we
rewrite (Φ/η)1/h as


Φ
η

1/h
=

Φ
.
η 1/h Φ1−1/h(8)Fig. 5. A comparison of the effective viscosity
for three plasticity models. While the effective
viscosity of a viscoplastic material is independent of stress, the effective viscosity of a shear
thinning Herschel-Bulkley material shrinks withincreasing stress. We also note that the effective
viscosity increases as a function of stress (see Appendix A) in the model of Bargteil et al. [2007].

This expression shows that the denominator η 1/h Φ1−1/h is theeffective viscosity coefficient from the vantage point of the standard
viscoplastic model (6); thus, a shear thinning material such as ours
has an effective viscosity that decreases as stress increases (seeFigure 5). When a larger stress is imposed, the foam flows moreeasily.
The Herschel-Bulkley model captures a wide range of viscoplastic
behaviors. While we use h < 1 for foams, the model can also
capture shear thickening behaviors as exemplified by materials suchas oobleck, a mixture of cornstarch and water (for a related treatment
via finite elements, see earlier work by Bargteil et al. [2007]).
Each of the increasingly sophisticated viscoplastic models wehave discussed generalizes its predecessor. Setting h = 1 recovers
the viscoplastic (Bingham) model; taking the limit η → 0 recovers
perfect plasticity.
Temporal Evolution of the Left Cauchy-Green Tensor. The

temporal evolution of the left Cauchy-Green tensor be typicallyinvolves an elastic prediction that applies the elastic deformation
followed by a plastic correction. Substituting our flow rule (5) into(4), we obtainḃe = Lbe + be LT −2
Tr[be ]γ(s)ŝ .
3(9)

We treat Eq. (9) through operator splitting, which first integrates theelastic predictionḃe = Lbe + be LT(10)followed by a plastic correction2ḃe = − Tr[be ]γ(s)ŝ .
3

5.2

(11)

Herschel-Bulkley Temporal Discretization

For numerical simulations, we must temporally discretize (10) and(11). Because foam and more generally Herschel-Bulkley materials undergo large shears, our discretization departs from that of
Stomakhin et al. [2013]. We extend the return mapping method, described by Simo and Hughes [1998] in the context of viscoplasticity,
to accommodate the Herschel-Bulkley power law (7). Within thiseframework it is sufficient to track b , but to implement our tearing
model, we additionally track F.
Elastic Prediction. At time step n + 1, we compute the incremen-tal deformation gradientfn+1 = (I + ∆t∇v n )(12)

due to all factors (e.g., elasticity, gravity) other than plasticity; v nis the particle velocity, and ∆t is the time step. We refer to the workof Stomakhin et al. [2013] for details on the computation of ∇v n .
Next we discretize (10) and predict a change only in the elasticstrain (see Appendix C), yielding an update to the volume preserving

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Continuum Foam: A Material Point Method for Shear-Dependent Flowseelastic strain bn+1 of
e T

e,pre

(13)

bn+1 = fn+1 bn fn+1 .

Likewise, we update Fn+1 via Fn+1 = fn+1 Fn .
Since s is the deviatoric part of the Kirchhoff stress tensor τ ,ewe have s = µ dev[b ], and the corresponding predicted stress
e,pre
pre
is sn+1 = µ dev[bn+1 ]. If the predicted stress’ magnitude spren+1
satisfies the yield condition (3), we accept the prediction.
Plastic Correction. If the yield condition is violated, we correct

the strain by discretizing (11) using fully implicit backward Euler.
Given the condition that volume does not change during plastic
correction, the discretization of (11) is given by
2
e
e,pre
e
bn+1 − bn+1 = − ∆t Tr[bn+1 ]γ(sn+1 )ŝn+1 .
3•

7

Algorithm 1 Updating the left Cauchy-Green tensor (9)Compute the change in the deformation gradient fn+1 using (12)Predict the volume-preserving part of the elastic strain using (13)Compute the predicted Kirchhoff stress using (2)if the yield condition is not violated (Φ ≤ 0) thenee,prebn+1 ← bn+1 and returnelse
Solve the nonlinear equation (16) for sn+1eUpdate bn+1 using equation (18)eeRenormalize bn+1 to ensure that det(bn+1 ) = 1
end if(14)

As shown in Appendix D, we transform this equation into a form
that requires the solution of a single scalar formulasn+1 − spre
n+1 = −2µ̃∆tγ(sn+1 ) ,

(15)

e,pre

where µ̃ = 13 Tr[bn+1 ]µ, or equivalently
η

1/h


sn+1 − spren+1 + 2µ̃∆tr
sn+1 −

2
σY
3

!1/h
= 0 . (16)This equation is well-defined as η → 0. When η = 0 or h = 1, we
obtain sn+1 explicitly as
! 
r

.
2
η
pre
pre
sn+1 = sn+1 − sn+1 −σY. (17)
1+3
2µ̃∆t
Otherwise, we solve Eq. (16) numerically. Denote the left hand sideof Eq. (16) as a function g(sn+1 ). The yield condition is violated
when we perform a plastic correction which implies that spre
n+1 >
q
q
pre
2
2
σ , hence g( 3 σY ) < 0 and g(sn+1 ) > 0. Additionally,
3 Yqg(sn+1 ) is monotonic in the interval 23 σY ≤ sn+1 ≤ spren+1 . We
thus solve for the unique root sn+1 via bisection.
Finally, we complete the correction by updating the Kirchhoff
stress tensor via sn+1 = sn+1 ŝpren+1 and by updating the volumetricleft Cauchy-Green strain tensor via
ebn+1 =

11e,presn+1 + Tr[bn+1 ]I .
µ
3

(18)

e

Note that when we update bn+1 in this way (i.e., by shrinking thedeviatoric part while keeping the trace part the same), the determinant of the result will always have a value equal to or slightly largerethan 1. We thus renormalize bn+1 after the update.
We summarize this numerical update in Algorithm 1.

6.

AVOID A VOID

Particle Resampling Under Large Deformation. In typicalfoams, the bulk modulus dramatically exceeds the shear modulus:
according to Weaire and Hutzler [2001], the shear modulus is on theorder of 101 Pa, while the bulk modulus is on the order of 105 Pa.
As a result, foams can undergo substantial shearing deformations.
When material points are numerically advected through large
deformations, their distribution can become highly nonuniform. For
example, when a sample of foam is compressed along one axis andFig. 6. Efficacy of our particle resampling scheme: (Left) A sheared
blob of foam develops voids if additional material points are not introduced
to the simulation. (Right) Our material point insertion technique preventsthe formation of voids.

stretched along a perpendicular axis, the uniformity of the material
point distribution can rapidly degrade, producing a sparse sampling
along the stretched axis. In a method for which inter-particle distancegoverns material connectivity, such as FLIP or MPM, this sparsity
results in the formation of numerically induced, nonphysical voids(Figure 6). In order to maintain a uniform point distribution, wedevelop a resampling approach for the material point method. As
described below, we employ Poisson disk sampling to insert new
points (§6.1), and we detect and merge points that are too close
to one another (§6.2). We further exploit an estimated distancefunction to restrict particle resampling to interior regions, effectively
distinguishing true exterior regions from artificial interior voids.

6.1

Inserting Points

We insert new material points to prevent the formation of numerically induced voids. We first compute an approximate narrowband
signed distance function (SDF) in the vicinity of the foam-air interface, enabling us to classify points as foam or air. We then identify
undersampled regions of the interior and insert new points to increase the sampling density.
Estimating the Interior Region. We compute the signed

distance function [Sethian 1999] on a uniform grid (distinct from the MPM background grid) of cell size r =1(cell width of the background grid). For the purpose of this op2
eration, we consider each material point to be a sphere of radius r,
approximating the foam volume as a union of spheres. By convention, the signed distance of a grid point to the surface of a sphere is

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

8•

Y. Yue et al.

Algorithm 2 Compute Signed Distance
1: S ← ∅. Initialize interface as empty sample set
2: for each grid point p do3:
p.d ← ∞
4: end for
5: Register all spheres into a uniform grid for range queries
6: for each grid point p do
7:
S 0 ← all spheres satisfying dist(p, center) < 2r
8:p.d ← min distance over all spheres ∈ S 0 , minus r
9: end for
10: for each grid edge e do
11:
if sgn(e.p1.d) 6= sgn(e.p2.d) then
. Zero crossing
12:
S.append zeroCrossing(e.p1, e.p2)
13:
end if
14: end for
15: for each grid point p do
16:
p.d ← sgn(p.d) · ∞17: end for
18: for each zeroCrossing C in S do
19:
for each grid point p satisfying dist(p, C) < 4r do
20:
if dist(p, C) < |p.d| then21:
p.d ← sgn(p.d) · dist(p, C)
22:
end if23:
end for
24: end for
positive when outside the sphere, zero on the surface, and negative
when inside the sphere. Algorithm 2 consists of three steps: first, wecategorize each grid point as interior or exterior based on its sign.
We then find the interpolated zero crossing points along the gridedges, and treat these as samples of the interface. Finally, we use
these interface samples to compute and assign improved estimatesof the distance magnitude for all nearby grid points (within 4r).
While the SDF magnitude is only meaningful in the vicinity of theboundary, its sign is meaningful everywhere. A point with negative
distance is without a doubt inside the material.
Octree-Based Dart Throwing. Poisson disk sampling [Cook1986] provides a simple technique to produce uniform sample distributions, while avoiding the insertion of samples too close to existing
points. Our method is the first application of Poisson disk sampling to MPM simulation, to our knowledge. We apply the octree
method proposed by Ebeida et al. [2012] in 3D. Whereas Ebeida
et al. randomly draw a sample from the active domain in a probability proportional to the cell volume, we instead draw one sample
uniformly from each active octree cell; refer to Algorithm 3 forpseudocode.
Put simply, we take a coarse to fine sampling approach. We begin
with one, active, coarse cell, and we repeat the following until noactive cells remain:
For each active cell, we generate a position sample at random.
We retain the sample if it lies inside of the foam, as determined bythresholding the SDF, and if it lies outside of all material spheres,
determined by using a modified radius of αr, where α < 1. This
ensures that material points can slightly overlap, thus preventing
the formationof numerically induced voids. We employ a value of√1
α ≈ 23 + 100
for all simulations. We choose α to be slightly larger
than half the diagonal length of a unit cube. With this choice, if
the particles remain within their initial distribution (8 particles percell), we do not insert additional particles; if the distance between
the particles grows larger than the initial setting, however, we insert
new particles. To avoid inserting particles too close to the interface,

Algorithm 3 Create Particles with Poisson Sampling
1: Ω ← {set of all material points}
2: A ← {set with one active cell that covers entire domain}
3: while A 6= ∅ do
. While there are active cells
4:for c ∈ A do
5:
Uniformly sample a point x in cell c6:
if x is not covered by any sphere of radius αr then7:
if dist(x) < −2.2 · r then. Sample the interior
8:
Accept x as a new sample
9:
end if10:
end if11:
Append all children of c to A
12:
Remove c from A
13:
end for14:
for c ∈ A do
15:
if c is covered by any sphere at p ∈ Ω of radius αr then16:
Remove c from A
17:
end if18:
if c ∩ Ω = ∅. Cell is void of material points. Cell is small
19:
or width(c) < √2r3 then20:
if (6 ∃ corner x of c) dist(x) ≤ −αr then
21:
Remove c from A
22:
end if23:
end if24:
end for
25: end while
which would induce popping artifacts, we only insert points if theinterpolated signed distance value is below a threshold of −2.2r. We
empirically observe that thresholds larger than this value introduce
noticeable popping artifacts. The threshold should be as large as
possible to maximize the resampling region, however. We thendeactivate the cell and recursively activate its eight children.
Before terminating the current iteration, we deactivate cells that
we deem unnecessary. We deactivate cells that are fully covered by aset of material spheres, as well as cells that are sufficiently far fromthe interior and that are too small or empty.
Determining Physical Values. We need to assign physical val-ues to newly inserted material points. For conserved quantities (i.e.,mass m and volume V ), we redistribute values from the surrounding
points to the new point. For each introduced point, we first search
for its neighboring points (e.g., within a distance r). Given N neighboring points with mass mq and volume Vq , where q = 1, ...,PN , weset the mass and volume of the added point to m = N1+1 q mq
P
and V = N1+1 q Vq , respectively. We set the masses and volumes of the neighboring points to m0q = mq − N1+1 mq and
Vq0 = Vq − N1+1 Vq , respectively.
e

For field values that are not inherently conserved (v, F, b ), we
employ a mass-weighted interpolation. First, we compute a fieldvalue qi at each grid point i using all of the material points thatexisted before resampling asP
p wip mp qp
qi = P,(19)
p wip mpwhere p is the index of the material point, wip is the shape function,
mp is the mass of the point, and qp is the field value at the point.
Next, we compute the field value q̃p for the new point p asXq̃p =wip qi .
(20)

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

i

Continuum Foam: A Material Point Method for Shear-Dependent Flows•9

allows the foam to return to a state in which it once again resists
stretching.

7.1

Fig. 7. The importance of an explicit tearing model: (Left) Without
explicit tearing, a thin, non-physical tendril forms. (Right) Our explicit
tearing model permits the thin tendril to break.
ee

To ensure that det(b ) = 1, we renormalize b . To avoid avolume change in F, which would correspond to change in its determinant, after interpolating both F and its determinant J, we rescale
F to ensure that its determinant is equal to J. Specifically, given theinterpolated value J and the initial interpolated value of F denotedF̃, we set F = J 1/3 (det(F̃)−1/3 F̃).

6.2

Merging Points

To avoid an overly dense sampling, we merge points that are too
close to one another. Like Ando et al. [2012], for each point wecompute its closest neighbor, and if the separating distance is below
a threshold (0.03r in our implementation), we add the pair to a list.
We then greedily find a set of disjoint point pairs, and merge eachpair into a single new point.
To assign values to a new merged point, we simply sum the valuesof the conserved quantities, and compute a mass-weighted averagefor the remaining quantities. As in the preceding section, we againerescale F for consistency with J and we renormalize b .

7.

TEARING

When a specimen of foam is increasingly stretched, it will grow thinand eventually tear. Standard MPM approaches lack an appropriate
explicit tearing model, however: as long as particles remain close to
one another relative to the grid spacing, they will exchange forces
through the background grid and behave as a continuous body. As aresult, foam simulated with MPM tends to become excessively thin
before finally separating, and in the process forms long, unrealistic
threads as shown in Figure 7.
To address this shortcoming, we propose to explicitly model tearing in our foam simulations. Physically, tearing occurs when theconnections between incident bubbles are broken, yielding an effectanalogous to particle dislocation in crystal structures. We first identify regions that are in a “weakened” or torn state. In these regions
we directly modify the computed forces and the accumulation ofstrain so that the material no longer resists separation, and so thatthe observed geometric stretching (i.e., the separation) does notaccumulate as additional physical strain. Finally, after tearing has
occurred, we gradually adjust the accumulated plasticity to model
the recovery of bubble neighbor connectivity in these regions. This

Detecting Weak Particles

We first iterate over all particles and label those whose accumulated
plasticity exceeds a threshold. Given a particle, its accumulated
plasticity is bp = Fp (Fp )T , where Fp is the plastic deformation
gradient. Similar to the yield condition, a particle should be labeled as weak if k dev[bp ]kF > σT , where the tear yield σT is aphenomenological material threshold.
To perform this check we need to compute k dev[bp ]kF . For
each particle, since we explicitly track the deformation gradienteF, its determinant J, and the normalized elastic strain b , we canepecompute k dev[b ]kF as follows. First, we recover b = J 2/3 b .
Next, denoting the right Cauchy-Green plastic deformation tensor
by Cp = (Fp )T Fp , we compute it as Cp = FT (be )−1 F. Finally,
the Frobenius norm of the left and right Cauchy-Green tensors are
the same, yielding
k dev[bp ]kF = k dev[Cp ]kF .

We compare this scalar quantity to σT to identify weakened particles.

7.2

Modeling the Weakening EffectOnce a region of foam is weakened, it should no longer resist separation and it should freely tear; weakened foam, however, should
continue to resist compression. To model this effect, we modify
the stress σ (computed in §5). First, we compute the eigendecomposition σ = VΣVT , where V is an orthonormal matrix, andΣ = diag(s1 , s2 , s3 ) is a diagonal matrix. We then set each eigenvalue si to min(0, si ), and reassemble σ. This modification sets
the stress corresponding to expansion to 0, while permitting thematerial to resist compression. (This is loosely similar to the workof Wang et al. [2010], who clamp the singular values of strain to
perform isotropic strain limiting.)
In addition, when we update the deformation gradient, we want toprevent regions that are in a weakened state from accumulating additional strain since we view them as torn. We therefore modify theincremental deformation gradient f in Eq. (12) to prevent the shear
from increasing and to prevent volumetric expansion, as follows.
To detect if material is accumulating shear, we consider the shear
before and after the proposed deformation. (Note that because thepreexisting shear and the incremental shear may have different directions, the incremental shear alone does not provide the necessary information.) We compute the eigendecompositions of the normalizedeleft Cauchy-Green tensor for the previous time step, bn , and of thee,preestimated left Cauchy-Green tensor for the current time step, bn+1
(from Eq. (13)), and obtain their maximum eigenvalues λn and
corr
pre
λpre
n+1 , respectively. If λn+1 ≤ λn , we let fn+1 = fn+1 . Otherwise,
the material is accumulating shear, so we compute Jf = det(fn+1 ),−1/3
we find the singular value decomposition fn+1 = Jf
fn+1 , andcorr
we set its singular values to 1 to obtain a corrected version fn+1 . We
then rescale to obtain the correct incremental deformation gradient
1/3 corr
fcorr
n+1 = Jf fn+1 . This effectively disables shear deformations
without affecting rotations.
Next, we disable the accumulation of expansion. We computecorr
Jn = det(Fn ) and Jfcorr = det(fcorr> 1,n+1 ). If Jn > 1 and Jf
this indicates that the material has a volume greater than its rest stateand is continuing to expand. In this case we simply normalize fcorr
n+1
to remove the expansion.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

10•

Y. Yue et al.

Finally, we proceed to use fcorr
n+1 in place of fn+1 in the steps that
follow Eq. (12).

7.3

Modeling Connectivity Recovery

After foam has torn, we do not want the material to stretch without
resistance indefinitely. Instead, we model the bubbles comprising
the foam as establishing connections with neighbors in the new
configuration, with the strength of these connections gradually recovering over a finite time scale. We consider the new state to bea bubble rearrangement without dislocation and therefore without
plasticity effects. In other words, we allow bp to gradually reduceto the identity tensor over the time scale of the rearrangement.
We therefore need a method to update F that incorporates thischange in the plasticity. Effectively, this requires manipulating Fp ,yet our particles do not explicitly carry this quantity. We will instead
compute Fp from Cp , which we can determine for each particle asin §7.1: Cp = Fp T Fp = FT (be )−1 F. While Cp is symmetric, Fp
is in general not symmetric, and we can thus only estimate Fp up
to some unknown rotation, i.e., Fp = RB, where R is an arbitrary
orthonormal matrix, and B is a symmetric matrix. Fortunately, thisextra rotational freedom will not affect our update of the plasticity:
we will only need to modify B.
To construct the above factorization, we first compute the eigendecomposition Cp √= USUT . We then obtain B as the principal rootpof C via B = U SUT . Next, we let X = FB−1 , giving us thedesired factorization of F into F = XB.
√
To update the plasticity, we would like to replace B = U√ SUT
with a “relaxed” version. We relax B by replacing S =diag(λ1 , λ2 , λ3 ), where λ1 , λ2 , λ3 > 0, with a new diagonal matrix P = diag(λ̃1 , λ̃2 , λ̃3 ). We need to ensure that det(Fp ) = 1 after this update, which in turn requires that det(P) = 1. Additionally,√
regarding P as a function of time, we require that P(0) = S andP(∞) = I to allow the foam to steadily transition from the current
plasticity level to zero plasticity. To model this behavior, we first emphasize the fact that det(Fp ) = λ1 λ2 λ3 = 1 and, that for any real
α αvalue α, λα1 λ2 λ3 = 1. If we regard α as a function of time and ifα(t)we let λ̃i (t) = λi , then det(P(t)) = 1 is automatically satisfied.
For α(t), we want to select a function that ensures that α(0)√ =1and that α(∞) = 0. P(t) then further satisfies P(0) = S andP(∞) = I. Specifically, we choose α(t) = exp(−t/ηp ), where
we call ηp the plasticity relaxation coefficient. Combining theseexp(−t/ηp )ingredients, we obtain λ̃i (t) = λi
, which is the solution to˙
the differential equation λ̃i (t) = − η1p λ̃i (t)log λ̃i (t) with the initial√ exp(−t/ηp )condition λ̃i (0) = λi . In matrix form this is P(t) = S.
00We compute the updated F as F = XB = XUPUT , where√ exp(−∆t/ηp )P= S.
Without our plastic recovery model, tearing regions are permanently weakened, resulting in a “brittle” effect. Setting ηp to a large
value gives brittle-type fracture behavior, whereas setting ηp to asmall value gives a more “sticky” tearing behavior. Refer to Figure 8
for a comparison of various ηp values.

7.4

Subgrid Geometry Removal

When a piece of foam grows excessively thin, the underlying discrete
representation could, in the worst case, be a single particle wide.
The physics governing these severely slender slivers can be strongly
under-resolved when computed on the background grid, giving riseto artificially stiff thread-like artifacts as in Figure 8 (center). While

Fig. 8. “Brittle”-type fracture (left, ηp = 100) v.s. sticky tearing (center,ηp = 0.35). Subgrid geometry removal erases geometry whose physics is
not resolved by the MPM background grid (right, ηp = 0.35). The mass
reduction due to subgrid geometry removal is less than 0.28% over the entire
simulation.

Table I. Material Parameters

Material
shaving creams’more interiors’more exteriorpieoobleck
viscoplastic

density[kg / m 3 ]
77.7
50.0
50.0
275.0
1,000.0
1,000.0

κ
[kPa]
109.0
109.0
109.0
109.0
109.0
109.0

µ
[kPa]
0.29
0.08
50.00
1.60
11.20
11.20

σY
[Pa]
31.9
10.0
1,000.0
120.0
0.1
0.1

η

h

σT

ηp

27.2
16.0
0.1
5.0
10.0
10.0

0.22
0.43
1.00
0.27
2.80
1.00

217.5
15.0
0.3
10.0
1.0
1.0

0.35
0.25
0.50
0.30
0.30
0.30

Parameters for all materials, including shaving cream (Figure 1), the interior and exteriorof a toasted s’more (Figure 11), a pie (Figure 12), a shear thickening oobleck material(Figures 13, and 14), and a viscoplastic material (Figure 13) for comparison with theoobleck.

a fully adaptive MPM treatment could ameliorate these issues, wedefer the exploration of this topic to future work.
Instead, we present a simple method to remove such underresolved particles by examining the anisotropy of the local region ina manner similar to Yu and Turk [2013]. For each weakened particlep, we identify all neighbors N (p) within a radius of 2h, where h is
the grid spacing. We then compute a mass-weighted mean position
pp and a mass-weighted covariance matrix Vp for each weakened
particle p asPq∈N (p) w(xq − xp )mq xq
pp = P(21)q∈N (p) w(xq − xp )mq
and
2 2
T
q∈N (p) w(xq − xp ) mq (xq − pp )(xq − pp )P22q∈N (p) w(xq − xp ) mp

PVp =, (22)

where w is the weighting function employed by Stomakhin etal. [2013]. We compute the eigendecomposition of Vp and consider the ratio between the largest and the smallest eigenvalues,
where a smaller ratio indicates a more anisotropic distribution ofneighboring particles. If this ratio is beneath a threshold (10−5 in our
tests) or the largest eigenvalue is too small (smaller than 10−5 ), weremove the weakened particle from the simulation. Figure 8 (right)
illustrates the effect of subgrid geometry removal.

8.

RESULTS

We now apply our model to simulations of foam and foam-like materials, as well as other Herschel-Bulkley materials, such as “oobleck.”
We list the material parameters from our simulations in Table I.
The physical parameters for real foam depend on a number of factors,
including freshness, subtle differences in formulation, temperature,etc., and hence can vary from product to product and from time totime. While the material density and the bulk modulus are rather

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Continuum Foam: A Material Point Method for Shear-Dependent Flowseasy to set, there are a number of key parameters that dictate theresulting motion, including the shear modulus µ, the yield stress σY ,the Herschel-Bulkley power h, the viscosity η, the tearing yield σT ,and the plastic relaxation ηp . We set the physical parameters (shear
modulus µ, yield stress σY , Herschel-Bulkley power h, viscosityη) of shaving cream based on rheological measurements from theliterature [Durian et al. 1990; Ovarlez et al. 2010]. In order to achieve
more compelling motion or to closely match the dynamics of aparticular physical sample, we adjusted some of these parameters
by at most 40% from the reported values, which is reasonable given
the likelihood of differences in our physical samples from those
employed in the rheology literature.
For materials whose true physical parameters are difficult to obtain, we were able to approximate any given desired behavior bytuning the six parameters according to the following procedure. We
begin by tuning the parameters µ and σY , which determine thebehavior of the elastic regime, before tuning h and η to adjust the
behavior of the plastic flow. Finally, we adjust σT and ηp to control
the tearing. The accompanying video provides an example of thisparameter tuning process.
Simulation statistics are given in Tables II and III. We timed all
examples on an 8-core Intel Xeon E5-2680 v2 2.8GHz processor.
Particle resampling is performed every 50 time steps. We show theamortized costs per time step in Fig. 9. At present we employ an
explicit integrator for simplicity, but we note that a semi-implicit
time stepping method is straightforward to implement by taking thechain rule with respect to be as derived in Stomakhin et al. [2013].
Finally, we employ a level-set based skinning method to extract
surfaces from the particles [Bhattacharya et al. 2011]. Please refer
to our video for the animation results.

Shaving Foam

11

Rasterize particle data to grid
17.3%

Particle resampling
35.6%

Tear detection

Particle positions
Particle-based collisions
Subgrid geometry removal
1.5%
Plastic recovery
13.4%

•

Plastic correction
1.7%

2.1%
Computeand limit grid forces
11.5%
Grid velocities
Grid-based collisions
1.2%
Update and limit
particle deformation data
9.8%
Particle velocities
6.1%

Particle resampling
12.6%
Particle-based collisions
Particle positions
1.4%

Pie Throw
Rasterize particle
data to grid
25.8%

Plastic recovery
15.0%
Tear detection
2.9%

Plastic correction
1.9%
Particle velocities
8.4%
Update and limit
particle deformation data
13.8%

Compute and limit
grid forces
16.1%
Grid velocities
Grid-based collisions
2.2%

Fig. 9. Timing breakdowns for a typical step in two simulations.

Shaving Foam. To test the validity of our model, we attached asample of real shaving foam to the base of a platform.
We acquired an approximation of the shape of the foam using
a 3D scanner, and subsequently recorded the motion of the foam
when subject to an oscillatory motion. We then simulated the acquired geometry, and compared the motion of our simulations to
the recorded video. Initially, the foam oscillates and slowly begins
to stretch downwards due to accumulated plasticity. As the oscillations continue, the foam exhibits typical nonlinear shear thinningbehavior, and rapidly tears apart. When the foam finally collides
with the ground plane, it exhibits a characteristic jiggle. As shownin the supplemental video and in Figure 1, our simulations exhibit
the same characteristic behaviors of real foam noted above.
The maximum volume deviation of the skinned shaving creammodel was 2.7%. Note that the density estimation (and the volume
estimation) can be erroneous near the interface. We investigated thedensity deviation for the interior of the foam. The deviation of theTable II. Simulation Statistics
Exampleshaving cream
uniform s’more
toasted s’morepieviscoplastic sphere
oobleck sphereviscoplastic penguin
oobleck penguin
oobleck pachinko#points
54,695- 72,681
202,278- 203,053
981,554- 988,355
725,052- 726,853
519,171- 534,871
519,171- 529,365
83,340- 87,030
83,340- 85,716
0-1,615,550

gridresolution
84× 84× 84
88×135× 88
140× 85×140
225×225×225
157×157×157
157×157×157
209×209×209
209×209×209
389×662× 56

min
frame
57.5
111.9
160.2
245.7
303.0
303.0
84.3
85.0
324.6

dt
1. 0 × 10−5

1. 0 × 10−5
1. 0 × 10−5
1. 0 × 10−5
0. 5 × 10−5
0. 5 × 10−5
0. 5 × 10−5
0. 5 × 10−5
0. 5 × 10−5

subgridgeom. rem.
X-

Statistics for the shaving cream shake (Figure 1), the toasted s’more (Figure 11), thepie toss (Figure 12), the viscoplastic and oobleck spheres (Figures 13 top row), theviscoplastic and oobleck penguins (Figure 13 bottom row), and the oobleck penguin
pachinko game (Figure 14). A frame is 1/30 of a second.

Fig. 10. S’more: A marshmallow made from a uniform material is crushed
between two graham crackers and begins to ooze through the cracks. We
then lift the upper cracker, resulting in the formation of a characteristic stiff
peak.

average density is at most 0.89% from the initial material densityof 77.7. The minimum and maximum densities were 18.0 and 93.5,
respectively. 95.7% of the grid density values were within 5% of
the average density.
Pie Throw. In Figure 12, we throw a whipped cream pie at aface. When the pie collides with the face, the collisions induce large
internal stresses, and the whipped cream flows. The pie subsequently
Table III. Grid Statistics
Exampleshaving cream
uniform s’more
toasted s’morepieviscoplastic sphere
oobleck sphereviscoplastic penguin
oobleck penguin
oobleck pachinko

cellspacing [mm]
3.0
1.6
1.0
1.6
1.6
1.6
1.2
1.2
1.8

initial#points/cell
8
8
8
8
8
8
8
8
8

#non-empty
cells
467,788-490,813
37,755- 53,529
133,806-164,153
135,628-181,445
85,455- 97,751
85,743- 90,629
21,559- 25,295
21,167- 23,681
0-461,329

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

12•

Y. Yue et al.

Fig. 13. Side by side comparisons between viscoplastic materials (left object in each pair) and shear thickening (oobleck) materials (right object in each pair).
While the viscoplastic materials flow immediately after impact with the ground plane, the shear thickening materials rebound elastically after the fast initial
impact, and only flow after the second, lower-velocity impact.

Fig. 11. Multi-material s’more: A toasted marshmallow composed ofa crunchy exterior and a gooey interior is crushed between two graham
crackers. The stiff outer layer cracks, allowing the softer interior to flow out.

Fig. 12. Pie throw: A whipped cream pie is thrown at a face and exhibits
the characteristic clumping and tearing behavior of such a system.

accumulates plasticity, stretches, and tears. Afterwards, chunks ofwhipped cream adhere to the face.
S’more. In Figure 10, we simulate a s’more by squeezing a marsh-

mallow between two graham crackers and a block of chocolate. The
marshmallow is modeled as a “moderately warm” foam, that is, as asoft material (we do not simulate the heat transfer). As we compress
the marshmallow, it begins to deform and flow, eventually seeping
through the holes in the graham cracker. As we subsequently lift theupper cracker, the marshmallow forms a thin neck and eventuallytears, leaving a characteristic stiff peak.
In Figure 11 we simulate a more complex marshmallow. We
model the marshmallow using two layers of MPM particles; thisallows us to mimic a toasted marshmallow, where a stiff outer layer
approximates the crispy exterior, and a soft inner layer approximatesthe molten marshmallow core (again, we do not simulate the heattransfer). After a critical amount of compression, the outer layer
cracks, allowing the softer inner layer to flow outwards.
Oobleck. Inspired by Bargteil et al. [2007], we also explore theapplication of our method to shear thickening materials such as“oobleck,” a solution of cornstarch in water. In Figure 13, we show a

Fig. 14. Oobleck pachinko: Penguins composed of shear thickening
oobleck are dropped onto a set of pachinko pins. The penguins bounce
after initial high-stress collisions, but later experience substantial plastic flow
after low-stress collisions.

side by side comparison between viscoplastic and shear thickening
materials. In this comparison, we only modify the Herschel-Bulkley
power parameter h; other parameters are the same for both materials.
We drop both materials onto a ground plane; upon impact theyexperience a large stress. At large stresses, the flow rate of the shear
thickening material is low, and it behaves like an elastic material andbounces off the ground. In contrast, the flow rate of the viscoplastic
material is much higher and it flows immediately.
In Figure 14, we use the same oobleck material to simulate apachinko machine. Due to shear thickening, the rapidly falling penguins initially bounce off of the pachinko pins, with only an occasional wing or foot being ripped away. As the penguins slow down
due to additional collisions with the increasingly tightly spaced pins,
the resulting gentle collisions induce significantly lower stresses,
and the viscoplastic penguins flow and lose their shape.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Continuum Foam: A Material Point Method for Shear-Dependent Flows9.

CONCLUSIONS AND FUTURE WORK

We have presented a method to simulate dense foams and othercomplex materials that exhibit nonlinear, viscoplastic behaviors. In
particular, we adapted the highly flexible Herschel-Bulkley constitutive model to the similarly powerful MPM framework in order torealistically simulate such materials. To robustly treat large shearing
effects characteristic of dense foams, we developed a particle resampling technique for MPM to prevent the formation of nonphysicalvoids. Finally, we extended MPM to treat tearing effects and to avoid
the generation of artificial thin strands.
While we successfully simulated dense foams and validated ourmodel against experimental results, there remain several opportunities for future research. First, it would be interesting to consider
additional types of foams. Our method simulates dense foams as acontinuum, homogenizing over a vast number of small bubbles inorder to conduct simulations at macroscopic scales. For less dense
foams where some individual bubbles are distinguishable (e.g., soap
foams in a bathtub) hybrid continuum-discrete strategies for both
simulation and rendering may be required, in order to achieve theappearance of high density bubbles without extreme computationalcosts.
Many real foams also exhibit additional phenomena such as coarsening and drainage, which modify the observed bulk properties ofthe material. Although the time scales of these phenomena are oftenmuch longer than our simulated time scale, it would be interesting
to incorporate such time-dependent material changes.
Finally, since a systematic rheological understanding of manyfoam phenomena is still lacking, we hope that the simulation tools
we have developed may aid in studying these phenomena from anumerical perspective.

Acknowledgements
We thank Keenan Crane for providing his digital and real face for thepie throw example, as well as Henrique Teles Maia, Nora Wixom andMichelle Ming-Yen Lee for helping to prepare the final version ofthis paper and the supplementary video. We thank Intel for donating
computing hardware, The Foundry for donating MODO licenses,
and Adobe for donating Creative Cloud licenses.
This work was supported in part by the JSPS Postdoctoral Fellowships for Research Abroad, NSF (Grants IIS-13-19483, CMMI-1129917, CAREER-1453101), NSERC (Grant RGPIN-04360-2014),
Intel, The Walt Disney Company, Autodesk, Side Effects, NVIDIA,
Adobe, and The Foundry.
REFERENCESA DAMS , B., PAULY, M., K EISER , R., AND G UIBAS , L. J. 2007. Adaptively
Sampled Particle Fluids. ACM Trans. Graph. 26, 3 (jul), 48:1–48:7.
A NDO , R., T H ÜREY, N., AND T SURUNO , R. 2012. Preserving Fluid Sheetswith Adaptively Sampled Anisotropic Particles. IEEE Transactions on
Visualization and Computer Graphics 18, 8 (aug), 1202–1214.
BANERJEE , B. 2004. Material Point Method Simulations of Fragmenting
Cylinders. In Proc. of 17th ASCE Engineering Mechanics Conference.
BARGTEIL , A. W., W OJTAN , C., H ODGINS , J. K., AND T URK , G. 2007.
A Finite Element Method for Animating Large Viscoplastic Flow. ACM
Trans. Graph. 26, 3 (jul), 16:1–16:8.
B HATTACHARYA , H., G AO , Y., AND BARGTEIL , A. 2011. A Level-set
Method for Skinning Animated Particle Data. In Proc. of ACM SIGGRAPH/Eurographics Symposium on Computer Animation 2011. 17–24.
B INGHAM , E. C. 1922. Fluidity and Plasticity. McGraw-Hill, 219.

•

13

B RACKBILL , J. U. AND RUPPEL , H. M. 1986. FLIP: A Method for Adaptively Zoned, Particle-in-Cell Calculations of Fluid Flows in Two Dimensions. J. Comp. Phys. 65, 2, 314–343.
B RAKKE , K. A. 1992. The Surface Evolver. Experimental Mathematics 1, 2,
141–165.
B USARYEV, O., D EY, T. K., WANG , H., AND R EN , Z. 2012. Animating
Bubble Interactions in a Liquid Foam. ACM Trans. Graph. 31, 4 (jul),
63:1–63:8.
C LEARY, P. W., P YO , S. H., P RAKASH , M., AND KOO , B. K. 2007.
Bubbling and Frothing Liquids. ACM Trans. Graph. 26, 3 (jul), 97:1–97:6.
C OHEN -A DDAD , S., H ÖHLER , R., AND P ITOIS , O. 2013. Flow in Foams
and Flowing Foams. Annual Review of Fluid Mechanics 45, 241–267.
C OOK , R. L. 1986. Stochastic Sampling in Computer Graphics. ACM Trans.
Graph. 5, 1 (jan), 51–72.
C OTTET, G.-H. AND KOUMOUTSAKOS , P. D. 2000. Vortex Methods:
Theory and Practice. Cambridge University Press.
D URIAN , D. J., W EITZ , D. A., AND P INE , D. J. 1990. Dynamics andCoarsening in Three-Dimensional Foams. Journal of Physics: Condensed
Matter 2, Supplement, SA433.
Ď URIKOVI Č , R. 2001. Animation of Soap Bubble Dynamics, Cluster Formation and Collision. Computer Graphics Forum 20, 3 (sep), 67–76.
E BEIDA , M. S., M ITCHELL , S. A., PATNEY, A., DAVIDSON , A. A., ANDOWENS , J. D. 2012. A Simple Algorithm for Maximal Poisson-Disk
Sampling in High Dimensions. Computer Graphics Forum 31, 2 (may),
785–794.
E DWARDS , E. AND B RIDSON , R. 2012. A High-Order Accurate Particlein-Cell Method. Int. J. Numer. Meth. Engng. 90, 9 (jun), 1073–1088.
E NRIGHT, D., F EDKIW, R., F ERZIGER , J., AND M ITCHELL , I. 2002. A
Hybrid Particle Level Set Method for Improved Interface Capturing. J.
Comp. Phys. 183, 1 (nov), 83–116.
G ARRETT, P. 1993. Recent developments in the understanding of foam
generation and stability. Chemical Engineering Science 48, 2, 367 – 392.
G ERSZEWSKI , D., B HATTACHARYA , H., AND BARGTEIL , A. W. 2009.
A Point-Based Method for Animating Elastoplastic Solids. In Proc. of
ACM SIGGRAPH/Eurographics Symposium on Computer Animation 2009.
133–138.
G OKTEKIN , T. G., BARGTEIL , A. W., AND O’B RIEN , J. F. 2004. A
Method for Animating Viscoelastic Fluids. ACM Trans. Graph. 23, 3(aug), 463–468.
G REENWOOD , S. T. AND H OUSE , D. H. 2004. Better with Bubbles:
Enhancing the Visual Realism of Simulated Fluid. In Proc. of ACM
SIGGRAPH/Eurographics Symposium on Computer Animation 2004. 287–
296.
G UO , Y. J. AND NAIRN , J. A. 2006. Three-Dimensional Dynamic Fracture Analysis Using the Material Point Method. Computer Modeling in
Engineering & Sciences 1, 1, 11–25.
H ARDER , A. AND M ANGNALL , C. 2013. Bubbles and Foam in Partysaurus
Rex. In ACM SIGGRAPH 2013 Talks. 17:1–17:1.
H ERSCHEL , W. H. AND B ULKLEY, R. 1926. Konsistenzmessungenvon Gummi-Benzollösungen. Kolloid-Zeitschrift & Zeitschrift für Polymere 39, 4 (aug), 291–300.
H IEBER , S. E. AND KOUMOUTSAKOS , P. 2005. A Lagrangian Particle
Level Set Method. J. Comp. Phys. 210, 1 (nov), 342 – 367.
H ONG , J.-M. AND K IM , C.-H. 2005. Discontinuous Fluids. ACM Trans.
Graph. 24, 3 (jul), 915–920.
H ONG , J.-M., L EE , H.-Y., YOON , J.-C., AND K IM , C.-H. 2008. Bubbles
Alive. ACM Trans. Graph. 27, 3 (aug), 48:1–48:4.
I HMSEN , M., A KINCI , N., A KINCI , G., AND T ESCHNER , M. 2012. Unified
Spray, Foam and Air Bubbles for Particle-Based Fluids. The VisualComputer 28, 6-8 (jun), 669–677.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

14•

Y. Yue et al.

I RVING , G., T ERAN , J., AND F EDKIW, R. 2004. Invertible Finite Elements for Robust Simulation of Large Deformation. In Proc. of ACM
SIGGRAPH/Eurographics Symposium on Computer Animation 2004. 131–
140.
J ONES , B., WARD , S., JALLEPALLI , A., P ERENIA , J., AND BARGTEIL ,A. W. 2014. Deformation Embedding for Point-Based Elastoplastic
Simulation. ACM Trans. Graph. 33, 2 (mar), 21:1–21:9.
K EISER , R., A DAMS , B., G ASSER , D., BAZZI , P., D UTR É , P., AND G ROSS ,M. 2005. A Unified Lagrangian Approach to Solid-Fluid Animation. In
Proc. of the Second Eurographics/IEEE VGTC Conference on Point-Based
Graphics 2005. 125–133.
K HAREVYCH , L., M ULLEN , P., OWHADI , H., AND D ESBRUN , M. 2009.
Numerical Coarsening of Inhomogeneous Elastic Materials. ACM Trans.
Graph. 28, 3 (aug), 51:1–51:8.
K IM , B., L IU , Y., L LAMAS , I., J IAO , X., AND ROSSIGNAC , J. 2007. Simulation of Bubbles in Foam with the Volume Control Method. ACM Trans.
Graph. 26, 3 (jul), 98:1–98:10.
K IM , D., S ONG , O.-Y., AND KO , H.-S. 2010. A Practical Simulation of
Dispersed Bubble Flow. ACM Trans. Graph. 29, 4 (jul), 70:1–70:5.
KOEHLER , S. A., S TONE , H. A., B RENNER , M. P., AND E GGERS , J. 1998.
Dynamics of foam drainage. Phys. Rev. E 58, 2097–2106.
K ÜCK , H., VOGELGSANG , C., AND G REINER , G. 2002. Simulation and
Rendering of Liquid Foams. In Proc. of Graphics Interface 2002. 81–88.
L OSASSO , F., S HINAR , T., S ELLE , A., AND F EDKIW, R. 2006. Multiple
Interacting Liquids. ACM Trans. Graph. 25, 3 (jul), 812–819.
L OSASSO , F., TALTON , J. O., K WATRA , N., AND F EDKIW, R. 2008. TwoWay Coupled SPH and Particle Level Set Fluid Simulation. IEEE Transactions on Visualization and Computer Graphics 14, 4 (jul-aug), 797–804.
M ÜLLER , M., K EISER , R., N EALEN , A., PAULY, M., G ROSS , M., AND
A LEXA , M. 2004. Point Based Animation of Elastic, Plastic and Melting Objects. In Proc. of ACM SIGGRAPH/Eurographics Symposium on
Computer Animation 2004. 141–151.
M ÜLLER , M., S OLENTHALER , B., K EISER , R., AND G ROSS , M.
2005. Particle-Based Fluid-Fluid Interaction. In Proc. of ACM SIGGRAPH/Eurographics Symposium on Computer Animation 2005. 237–
244.
NARAIN , R., G OLAS , A., AND L IN , M. C. 2010. Free-Flowing GranularMaterials with Two-Way Solid Coupling. ACM Trans. Graph. 29, 6 (dec),
173:1–173:9.
N ESME , M., K RY, P. G., J E Ř ÁBKOV Á , L., AND FAURE , F. 2009. Preserving
Topology and Elasticity for Embedded Deformable Models. ACM Trans.
Graph. 28, 3 (aug), 52:1–52:9.
N IELSEN , M. B. AND Ø STERBY, O. 2013. A Two-Continua Approach
to Eulerian Simulation of Water Spray. ACM Trans. Graph. 32, 4 (jul),
67:1–67:10.
O’B RIEN , J. F., BARGTEIL , A. W., AND H ODGINS , J. K. 2002. GraphicalModeling and Animation of Ductile Fracture. ACM Trans. Graph. 21, 3(jul), 291–294.
OVARLEZ , G., K RISHAN , K., AND C OHEN -A DDAD , S. 2010. Investigation of shear banding in three-dimensional foams. EPL (Europhysics
Letters) 91, 6, 68005:1–68005:6.
PATKAR , S., A ANJANEYA , M., K ARPMAN , D., AND F EDKIW, R. 2013.
A Hybrid Lagrangian-Eulerian Formulation for Bubble Generation and
Dynamics. In Proc. of ACM SIGGRAPH/Eurographics Symposium on
Computer Animation 2013. 105–114.
S AYE , R. I. AND S ETHIAN , J. A. 2013. Multiscale Modeling of Membrane Rearrangement, Drainage, and Rupture in Evolving Foams. Science 340, 6133 (may), 720–724.
S CHECHTER , H. AND B RIDSON , R. 2012. Ghost SPH for Animating Water.
ACM Trans. Graph. 31, 4 (jul), 61:1–61:8.

S CHREYER , H. L., S ULSKY, D. L., AND Z HOU , S.-J. 2002. Modeling
Delamination as a Strong Discontinuity with the Material Point Method.
Computer Methods in Applied Mechanics and Engineering 191, 23–24(mar), 2483 – 2507.
S ETHIAN , J. A. 1999. Level Set Methods and Fast Marching Methods.
Cambridge University Press.
S IMO , J. C. 1988. A Framework for Finite Strain Elastoplasticity Based on
Maximum Plastic Dissipation and the Multiplicative Decomposition: Part
I. Continuum Formulation. Computer Methods in Applied Mechanics andEngineering 66, 2 (feb), 199–219.
S IMO , J. C. AND H UGHES , T. J. R. 1998. Computational Inelasticity.
Springer.
S OLENTHALER , B. AND PAJAROLA , R. 2008. Density Contrast SPH
Interfaces. In Proc. of ACM SIGGRAPH/Eurographics Symposium on
Computer Animation 2008. 211–218.
S OLENTHALER , B., S CHL ÄFLI , J., AND PAJAROLA , R. 2007. A Unified
Particle Model for Fluid-Solid Interactions. Computer Animation andVirtual Worlds 18, 1 (feb), 69–82.
S TEFFEN , M., K IRBY, R. M., AND B ERZINS , M. 2008. Analysis andReduction of Quadratire Errors in the Material Point Method (MPM). Int.
J. Numer. Meth. Engng. 76, 6 (nov), 922–948.
S TOMAKHIN , A., S CHROEDER , C., C HAI , L., T ERAN , J., AND S ELLE ,A. 2013. A Material Point Method for Snow Simulation. ACM Trans.
Graph. 32, 4 (jul), 102:1–102:9.
S TOMAKHIN , A., S CHROEDER , C., J IANG , C., C HAI , L., T ERAN , J.,
AND S ELLE , A. 2014. Augmented MPM for Phase-Change and Varied
Materials. ACM Trans. Graph. 33, 4 (jul), 138:1–138:11.
S ULSKY, D. AND S CHREYER , L. 2004. MPM Simulation of Dynamic
Material Failure with a Decohesion Constitutive Model. European Journalof Mechanics - A/Solids 23, 3 (may-jun), 423 – 445.
S ULSKY, D., Z HOU , S.-J., AND S CHREYER , H. L. 1995. Applicationof a Particle-in-Cell Method to Solid Mechanics. Computer Physics
Communications 87, 1-2 (may), 236–252.
T ERZOPOULOS , D. AND F LEISCHER , K. 1988. Modeling Inelastic Deformation: Viscoelasticity, Plasticity, Fracture. Computer Graphics 22, 4(aug), 269–278.
T H ÜREY, N., S ADLO , F., S CHIRM , S., M ÜLLER -F ISCHER , M., AND
G ROSS , M. 2007. Real-time Simulations of Bubbles and Foam within a
Shallow Water Framework. In Proc. of ACM SIGGRAPH/Eurographics
Symposium on Computer Animation 2007. 191–198.
WANG , H., O’B RIEN , J., AND R AMAMOORTHI , R. 2010. Multi-Resolution
Isotropic Strain Limiting. ACM Trans. Graph. 29, 6 (dec), 156:1–156:10.
W EAIRE , D. 2008. The Rheology of Foam. Current Opinion in Colloid &Interface Science 13, 3, 171–176.
W EAIRE , D. AND H UTZLER , S. 2001. The Physics of Foams. Oxford
University Press.
W ICKE , M., R ITCHIE , D., K LINGNER , B. M., B URKE , S., S HEWCHUK ,J. R., AND O’B RIEN , J. F. 2010. Dynamic Local Remeshing for Elastoplastic Simulation. ACM Trans. Graph. 29, 4 (jul), 49:1–49:11.
W OJTAN , C. AND T URK , G. 2008. Fast Viscoelastic Behavior with Thin
Features. ACM Trans. Graph. 27, 3 (aug), 47:1–47:8.
Y U , J. AND T URK , G. 2013. Reconstructing Surfaces of Particle-Based
Fluids Using Anisotropic Kernels. ACM Trans. Graph. 32, 1 (jan), 5:1–
5:12.
Z HENG , W., YONG , J.-H., AND PAUL , J.-C. 2006. Simulation of Bubbles.
In Proc. of ACM SIGGRAPH/Eurographics Symposium on Computer
Animation 2006. 325–333.
Z HOU , S. 1998. The Numerical Prediction of Material Failure Based on the
Material Point Method. Ph.D. thesis, The University of New Mexico.
Z HU , Y. AND B RIDSON , R. 2005. Animating Sand as a Fluid. ACM Trans.
Graph. 24, 3 (jul), 965–972.

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Continuum Foam: A Material Point Method for Shear-Dependent Flows•

15

DERIVATION OF EQ. (15)APPENDIXD.

A.

Taking the trace of both sides of Eq. (14) and noting that Tr[ŝn+1 ] =0, we find thatTHE EFFECTIVE VISCOSITY IN THE MODEL
OF BARGTEIL ET AL.

e

Bargteil et al. [2007] proposed the function

kPkF − Py
γ = min ν
,1 ,
kPkF

e,pre

sn+1 − spre
n+1 = −2µ̃∆tγ(sn+1 )ŝn+1 .
Equivalently, we havepresn+1 ŝn+1 − spren+1 ŝn+1 = −2µ̃∆tγ(sn+1 )ŝn+1 ,From (1), we have
∂W eT
F =
∂Fe

e

∂Wv (J)∂Ws (b )
+
∂Fe
∂Fe

!
FeT .

(23)

Next, using the chain rule, we obtain∂Wv (J)∂Wv (J) ∂J
1
1 ∂ det(Fe )== κ(J − )
e
e
∂F
∂J ∂F
2
J
∂Fe
1
1
1
= κ(J − ) det(Fe )Fe−T = κ(J 2 − 1)Fe−T .
2
J
2
(24)
Likewise,
e

e

e

e

∂Ws (b ) ∂ Tr[b ] ∂b
∂Ws (b )
=
:
,
e
e
∂Fe
∂Fe
∂ Tr[b ] ∂b

(25)

where we haveee∂Ws (b ) ∂ Tr[b ]e∂ Tr[b ]

∂b

e

=

1
µI
2

(26)

and
e

∂be ij∂b ij∂J −2/3= J −2/3 e + be ij
e
∂F uv
∂F uv
∂Fe uv= J −2/3 (Fe iv δju + Fe jv δiu ) −2 −2/3 e e−TJb ij F uv .
3
(27)

Substituting (26) and (27) into (25), we obtaine∂Ws (b )
1
−2/3
e
e e−T
=
µJ
F
−
Tr[b
]F
.
∂Fe
3

(28)

Finally, substituting (24) and (28) into (23), we obtain (2).

C.

e,pre
e
e
1
Tr[bn+1 ]I. Likewise, bn+1 = µ1 sn+1 + 31 Tr[bn+1 ]I = µ1 sn+1 +
3
e,pre
1
Tr[bn+1 ]I. Thus, Eq. (14) becomes
3

DERIVATION OF EQ.(2)

τ =

e,pre

1 pre
Recalling that spre
n+1 = µ dev[bn+1 ], we have bn+1 = µ sn+1 +where P gives the stress, Py gives the yield stress, and ν gives theflow rate. We have omitted work hardening and work softening forsimplicity. The effective viscosity is thus kPkF /ν, an increasing
function of the stress.

B.

e,pre

Tr[bn+1 ] − Tr[bn+1 ] = 0 .

DERIVATION OF EQ.(13)eee Tbe,pren+1 = bn + ∆t(Ln bn + bn Ln ) .

We can add a higher order term ∆t2 (Ln ben LTn ) to the right handside, which vanishes as ∆t → 0. We then obtain= (I + ∆tLn )ben (I + ∆tLTn ) = fn+1 ben fn+1 .
Normalizing both sides, we obtain Eq. (13).

E.

ALGORITHM REFERENCE

E.1

MPM Discretization

As in [Stomakhin et al. 2013], we use a regular grid as the background Eulerian mesh for the computation of stress-induced forces,
which requires the spatial discretization of ∇ · σ. We fix the grid
position and the grid spacing throughout the simulation to avoid
remeshing and to simplify computations.
We use cubic B-splines [Steffen et al. 2008; Stomakhin et al.
2013] for the grid basis functions Nih (xp ):1
1
1
Nih (xp ) = N ( (xp − ih))N ( (yp − jh))N ( (zp − kh))hhhHere i = (i, j, k) is the grid index, xp = (xp , yp , zp ) is the evaluation position (or the particle position), h is the grid spacing, and

1
2
3
2

0 ≤ |x| < 1
 2 |x| − x + 3 ,1
3
2
N (x) = − 6 |x| + x − 2|x| + 43 , 1 ≤ |x| < 2 .
0,otherwise
For compact notation, let wip = Nih (xp ) and ∇wip = ∇Nih (xp )as in [Stomakhin et al. 2013].
Physical quantities along the flow are tracked with Lagrangian particles, which are also useful for tracking history dependent quantities
required for the computation of plasticity. Such physical quantities
include the position x, the velocity v, the mass m, the total defore
mation gradient F, and the normalized strain b . A complete list of
the quantities stored for each particle can be found in §E.4.

E.2

Discretizing Eq. (10), we obtain

e
e
e T
2
e T
be,pren+1 = bn + ∆t(Ln bn + bn Ln ) + ∆t (Ln bn Ln )which reveals that ŝn+1 = ŝpren+1 . That is, s changes in magnitude
but not in direction. It thus suffices to solve the scalar equation (15).

Initial Set Up

As in [Stomakhin et al. 2013], we compute each particle’s volume
and density once during the firstPtime step. A particle’s density andvolume are estimated as ρp = i wip mi /h3 and Vp = mp /ρp ,respectively.

E.3

Per Step Computation

To update the physical quantities at time step n + 1 given those attime step n, we perform the following steps in order. For brevity,
we drop the subscript n + 1 and n when it is clear from the context.
ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

16•

Y. Yue et al.

1. Rasterize particle data to the P
grid. We transfer mass fromthe particles to the grid via mi =p wip mp , where mp is themass of particle p. Likewise, we transfer momentum via mi v i =P
v p is the velocity of particle p. This yields ap wip mp v p , where P
grid velocity of v i = p wip mp v p /mi .
2. Detect tearing regions. For each particle p, we computek dev[Cpp ]kF = k dev[FT (be )−1 F]kF . We label p as weak only
if k dev[Cpp ]kF > σT .
3. Compute grid forces. We compute the internal forces from thestress (∇ · σ) at each grid node i viaX
fi = −Vp Jp σp ∇wip ,pwhere Vp is the initial particle volume and Jp = det(Fp ) is thedeterminant of the particle total deformation gradient Fp . σp is
computed from the constitutive relation viai
1 hκ 2e
σp =(Jp − 1)I + µ dev[bp ] ,Jp 2ewhere bp is the normalized strain of particle p.
4. Limit forces. If a particle is labeled as weak in step 3, we modify its stress σp by first computing the eigendecomposition σp =Vp Σp VTp , where Vp is a unitary matrix and Σp = diag(s1 , s2 , s3 )is a diagonal matrix. Next, we replace each eigenvalue si bymin(0, si ), and reassemble σp .
5. Update grid velocities. We set the grid velocities v ? to

9. Update particle velocities. Following [Stomakhin et al. 2013],
we compute new particle velocities
as v p,n+1 = (1 − α)v PICp +P
αv FLIPp , where v PICp =wip v ?i and v FLIPp = v p,n +iP?
i wip (v i − v i ). We employ α = 0.95.

preplastic flow. We compute sn+1=e,preµk dev[bn+1 ]kF . If the yield condition is violated, i.e., if
Φ(spren+1 ) > 0, we compute the plastic flow by solving for sn+1

10.

Computein (16). When η = 0 or h = 1, an explicit solve is possible, withsn+1 given by (17). Otherwise, wehq
solve Eq. (16)inumerically via
bisection starting from the interval

2
σ , spre
n+1
3 Y

.

Next, we complete the correction by updating sn+1 = sn+1 ŝpren+1
and by updating the volumetric left Cauchy-Green tensor via (18).
eFinally, we renormalize bn+1 .
11. Compute plastic recovery. We first compute Cp=
FT (be )−1 F. Next,
the eigendecomposition Cp =√ we compute√
USUT to obtain Cp as U SUT . √Afterwards, we let λi = ( S)ii and compute λ̃i =exp(−∆t/ηp )λi
to obtain R = diag(λ̃i /λi ).
Finally, we update the deformation gradient to F0 = FURUT .
12. Remove subgrid geometry. For each weakened particle p,
we identify all neighbors N (p) within a radius of 2h, where h is
the grid spacing. We then compute a mass-weighted mean position
pp and a mass-weighted covariance matrix Vp for each weakened
particle p via (21) and (22), respectively. The weighting function is
given byv ?i = ∆tf i /mi + ∆tg .

1
1
1
w(xq − xp ) = N ( (xq − xp ))N ( (yq − yp ))N ( (zq − zp )) .
hhh

6. Resolve grid-based collisions. To model adhesive obstacles,
we set the relative velocity between a grid node and an obstacle to 0
if the grid node is inside an obstacle.
We model non-adhesive obstacles as in [Stomakhin et al. 2013]§8.

We compute the eigendecomposition of Vp and consider the ratio
between the largest and smallest eigenvalues. If this ratio is beneatha threshold (10−5 in our tests) or the largest eigenvalue is too small(smaller than 10−5 ), we remove the weakened particle from thesimulation.

7. Update particle deformation data. We compute the incremen-tal deformation gradientPvia fp,n+1 = (I+∆t∇v p,n ), where ∇v p,n
is given by ∇v p,n = i v i,n (∇wip )T .

T
e,pre
e
We compute the predicted strain via bp,n+1 = fp,n+1 bp,n fp,n+1 .

Likewise, we update Fp,n+1 to Fp,n+1 = fp,n+1 Fp,n .

8. Limit deformation gradient. Based on the results of step 7,

if a particle is weak and if the shear component is increasing, wedisable the shear deformation by modifying f. In addition, if thevolume continues to expand, we disable the expansion.
To detect accumulation, we compute the eigendecompositions ofee,prebp,n and bp,n+1 , and compare their maximum eigenvalues λn andprepre
λn+1 . If λn+1 ≤ λn , we set fcorr
n+1 = fn+1 . Otherwise, the material is accumulating shear, and we compute Jf = det(fn+1 ). We−1/3
compute the singular value decomposition of fn+1 = Jf
fn+1 ,
corr
setting its singular values to 1 to obtain a corrected version, fn+1 .
corr
1/3
We then rescale again to obtain fcorr
n+1 = Jf fn+1 .
To detect continued volume expansion, we compute Jn =corr
det(Fn ) and Jfcorr = det(fcorr> 1, then+1 ). If Jn > 1 and Jf
material was previously expanding and continues to expand. In thissituation, we disable this additional expansion by normalizing fcorr
n+1 .
e,pre
Finally, we update the predicted strain to bp,n+1 =ecorr Tfcorr
n+1 bp,n fn+1 .

13. Resolve particle-based collisions. To model adhesive ob-

stacles, we set the relative velocity between a particle and an obstacle
to 0 if the particle is inside an obstacle.
We model non-adhesive obstacles as in [Stomakhin et al. 2013]§8.
14. Update particle positions. We set the particle positions to

xp,n+1 = xp,n + ∆tv p,n+1 .

15. Avoid a Void. We fill voids for a sparse subset of all time

steps. We consider each material point to be a sphere of radiusr = 12 h, where h is the spacing of the MPM background grid.
We then employ Algorithm 2 to compute an approximate signed
distance function (SDF) on a uniform grid (distinct from the MPM
background grid) of cell size r. We fill potential void regions by
first reducing
the radius of each point to αr. We use a value of√1
α ≈ 23 + 100
in our simulations. We then employ octree-based
dart throwing (Algorithm 3) to populate interior regions (2.2r awayfrom the interface) with new particles.
We need to assign physical values to newly inserted material
points. For conserved quantities (i.e., mass m and volume V ), we
redistribute values from the surrounding points to the new point.
e
For field values that are not inherently conserved (v, F, b ), weeemploy mass-weighted interpolation. We renormalize b to ensure

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Continuum Foam: A Material Point Method for Shear-Dependent Flowsethat det(b ) = 1. After interpolating both F and its determinant J,
we rescale F to ensure that its determinant is equal to J.
Finally, we merge points that are too close to one another. To
assign values to a new merged point, we simply sum the values of
the conserved quantities, and compute a mass-weighted average forthe remaining quantities. We again rescale F for consistency with J,e
and renormalize b .

E.4

Internal Data Structure and Variables

MPM Background Grid. At each MPM grid node i, we store themass mi , the velocity v i , the updated velocity v ?i , and the internal
force f i . In addition, we allocate memory for the deformation
gradient Fi , for the determinant of the deformation gradient Ji , forethe strain bi , and for the sum of weights wi,sum (in order to compute
field values at newly created particles).
Material Points. Each material point p stores its mass mp , itsevolume Vp , its velocity v p , its normalized strain b , and its total
deformation gradient F. In addition, we store the binary flags weak
and valid. The flag weak indicates tearing, and is true if and only
if k dev[Cpp ]kF > σT . The flag valid indicates whether or not wewill remove a given particle.
Auxiliary Uniform Grid. We employ an auxiliary uniform gridto accelerate neighboring particle-particle queries. Note that otherneighbor queries do not require an auxiliary uniform grid. That is,neighboring grid node queries for particles can be directly computedby converting particle positions to grid node IDs, and neighboring
particle queries for grid nodes can be avoided by inverting thecomputation order: we scatter particle data onto the grid.
At each node of the uniform grid, we store a list of the neighborparticles.
SDF Grid. At each node of the SDF grid, we store the minimum

distance to the interface and a pointer to the closest zero-crossing
point.

E.5

Pseudocode

Algorithm 4 Simulate Foam
1: Compute Particle Volumes and Densities
2: time step ← 0
3: while time step < num steps do
4:
Rasterize Particle Data to Grid5:
Detect Tearing Regions
6:
Compute and Limit Grid Forces
7:
Update Grid Velocities
8:
Resolve Grid Based Collisions
9:
Update and Limit Deformation Gradient
10:
Update Particle Velocities
11:
Compute Plastic Flow
12:
Compute Plastic Recovery
13:
Remove Subgrid Geometry
14:
Resolve Particle Based Collisions
15:
Update Particle Positions
16:
if time step mod Avoid a Void Interval = 0 then
. Avoid a Void Interval = 50
17:
Avoid a Void18:
end if19:
time step ← time step + 1
20: end while•

17

Algorithm 5 Compute Particle Volumes and Densities
1: Rasterize Particle Data to Grid
2: for each particle p do3:
ρp ← 04:
Mp ← all MPM grid nodes satisfying wip > 0
5:for each MPM grid node i ∈ Mp do
6:
ρp ← ρp + wip mi. Accumulate mass
7:
end for8:ρp ← ρp /h3
9:Vp ← mp /ρp
10: end for

Algorithm 6 Rasterize Particle Data to Grid
1: for each MPM grid node i do2:mi ← 03:vi ← 04: end for
5: for each particle p do
6:
Mp ← all MPM grid nodes satisfying wip > 0
7:
for each MPM grid node i ∈ Mp do
8:mi ← mi + wip mp. Accumulate mass
9:v i ← v i + wip mp v p. Accumulate momentum10:
end for
11: end for
12: for each MPM grid node i do
13:v i ← v i /mi
14: end for

Algorithm 7 Detect Tearing Regions
1: for each particle p doe2:
Cpp ← det(Fp )−2/3 FTp (bp )−1 Fpp3:
if k dev[Cp ]kF > σT then4:
Label p as weak
5:
else
6:
Label p as not weak7:
end if
8: end for

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

18•

Y. Yue et al.

Algorithm 8 Compute and Limit Grid Forces
1: for each MPM grid node i do
2:fi ← 03: end for
4: for each particle p do
5:
Jp ← det(F
i
h p)
e

σp ← J1p κ2 (Jp2 − 1)I + µ dev[bp ]7:
if p is labeled as weak then8:Vp Σp VTp ← σp
. Eigendecompostion9:si ← (Σp )ii. Eigenvalues10:Σp ← diag(min(0, si )). Limit stress11:
σp ← Vp Σp VTp
12:
end if13:
Mp ← all MPM grid nodes i satisfying wip > 0
14:
for each MPM grid node i ∈ Mp do
15:
f i ← f i − Vp Jp σp ∇wip
. Accumulate force
16:
end for17: end for6:

Algorithm 9 Update Grid Velocities
1: for each MPM grid node i do
2:
v ?i ← ∆tf i /mi + ∆tg
3: end for

Algorithm 10 Resolve Grid Based Collisions
1: for each MPM grid node i do
2:
if i is inside an obstacle then3:
v ?i ← v obstacle4:
end if
5: end for

Algorithm 11 Update and Limit Deformation Gradient
1: for each particle p do
2:∇v p,n ← 0
3:
Mp ← all MPM grid nodes i satisfying wip > 0
4:for each MPM grid node i ∈ Mp do
5:∇v p,n ← ∇v p,n + v i,n (∇wip )T6:
end for7:
fp,n+1 ← (I + ∆t∇v p,n )
8:
9:
10:
11:
12:
13:

e,pre

e

T

bp,n+1 ← fp,n+1 bp,n fp,n+1
if p is not labeled as weak thenFp,n+1 ← fp,n+1 Fp,ncontinueend ife
λn ← largest eigenvalue of bp,n
e,pre
14:
λpre
n+1 ← largest eigenvalue of bp,n+1
corr
15:
fn+1 = fp,n+1
16:
if λpre. Shear is accumulatingn+1 > λn then17:
Jf ← det(fp,n+1 )−1/3
18:
fn+1 ← Jf
fp,n+1. SVD
19:USVT ← fn+1
20:
fn+1 ← UVT
T
e,pre
e
21:
bp,n+1 ← fn+1 bp,n fn+1
1/3
. Only retain the rotation
22:
fcorr
n+1 ← Jf fn+1
23:
end if24:
Jn = det(Fp,n )25:
Jfcorr = det(fcorr
n+1 )
26:
if Jn Jfcorr > 1 and Jfcorr > 1 then. Previously and continuing to expand
corr −1/3 corr
27:
fcorr
)
fn+1
. Discard expansion
n+1 ← (Jf28:
end if29:
Fp,n+1 ← fcorr
n+1 Fp,n
30: end for

Algorithm 12 Update Particle Velocities
1: for each particle p do
2:
v PICp ← 0
3:
v FLIPp ← v p
4:
Mp ← all MPM grid nodes i satisfying wip > 0
5:for each MPM grid node i ∈ Mp do
6:v PICp ← v PICp + wip v ?i
7:
v FLIPp ← v FLIPp + wip (v ?i − v i )8:end for9:v p ← (1 − α)v PICp + αv FLIPp. α = 0.95
10: end for

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

Continuum Foam: A Material Point Method for Shear-Dependent Flows

Algorithm 13 Compute Plastic Flow
1: for each particle p do
e,pre
2:
spre
n+1 ← µ dev[bn+1 ]
pre
pre
3:
sn+1 ← ksn+1 kF
q2
4:
if spreσ ≤ 0 then
n+1 −
3 Y
5:
continue6:
end ife,pre7:µ̃ ← 31 Tr[bn+1 ]µ8:if η = 0 or h = 1 then
q

9:

10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:

. Elastic regime .  . Direct solveη
1 + 2µ̃∆telse. Numerical solve via bisectionq. Initial search range[sm , sM ] ← [ 23 σY , spren+1 ]repeats ← (sm + sM )/2
q

 h1

1
2
g ← η h s − spre
+
2µ̃∆t
s
−
σ
Y
n+1
3
if g < 0 thensm ← s. Update search rangeelsesM ← s. Update search range
end ifE ← g/spre. Relative error
n+1
until |E| < .  = 10−6
end if
pre
pre
ŝpre
. Flow direction
n+1 ← sn+1 /sn+1
presn+1 ← sn+1 ŝn+1
ee,prebn+1 ← µ1 sn+1 + 13 Tr[bn+1 ]I
sn+1 ← spre
n+1 −

e

spre
n+1 −

e

e

26:
bn+1 ← det(bn+1 )−1/3 bn+1
27: end for

2
σ
3 Y

. Renormalize•

19

Algorithm 15 Remove Subgrid Geometry
1: for each particle p do
2:
if p is labeled as weak then3:
Np ← all particles q satisfying w(xp − xq ) > 0. Use uniform grid for query4:
pp ← 0. Initialize mean position5:
wp,sum ← 0. Initialize weight sum
6:for each particle q ∈ Np do
7:w ← w(xq − xp )mq. Compute weight8:
pp ← pp + wxq. Accumulate position9:
wp,sum ← wp,sum + w. Accumulate weight10:
end for11:pp ← pp /wp,sum. Compute mean position12:
Vp ← 0. Initialize covariance matrix
13:
wp,sum ← 0. Initialize weight sum
14:for each particle q ∈ Np do
15:w ← (w(xq − xp ))2 m2q. Compute weight16:
Vp ← Vp + w(xq − pp )(xq − pp )T. Accumulate covariance
17:
wp,sum ← wp,sum + w. Accumulate weight18:
end for19:
Vp ← Vp /wp,sum. Covariance matrix
20:
USUT ← Vp. Eigendecomposition
21:(λM , λm ) ← (maxi {Sii }, mini {Sii })22:
if λm /λM <  then.  = 10−5
. Neighbor distribution is highly anisotropic
23:
Label p as invalid
24:
end if25:
end if
26: end for
27: for each particle p do
28:
if p is invalid then29:
Remove p from the simulation
30:
end if
31: end for

Algorithm 16 Resolve Particle Based Collisions
1: for each particle p do
2:
if p is inside an obstacle then3:v p ← v obstacle4:
end if
5: end for

Algorithm 14 Compute Plastic Recovery
1: for each particle p doe2:
Cpp ← det(Fp )−2/3 FTp (bp )−1 Fp
3:
USUT ← Cpp. Eigendecomposition

Algorithm 17 Update Particle Positions
1: for each particle p do2:
xp ← xp + ∆tv p
3: end for(exp(−∆t/ηp )−1)/24:
λi ← Sii
5:
R ← diag(λi )6:
Fp ← Fp URUT
7: end for

Algorithm 18 Avoid a Void
1: Compute Signed Distance. Algorithm 2
2: Create Particles with Poisson Sampling. Algorithm 3
3: Compute New Conserved Values. Algorithm 19
4: Compute New Field Values. Algorithm 20
5: Merge Particles. Algorithm 21
ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.

20•

Y. Yue et al.

Algorithm 19 Compute New Conserved Values
1: for each newly sampled particle p do2:
Np ← all existing particles q satisfying |xq − xp | < r
3:
N ← |Np |. The number of neighbor particles
4:
mp,sum ← 0. Initialize total mass
5:
Vp,sum ← 0. Initialize total volume
6:
for each particle q ∈ Np do
7:mp,sum ← mp,sum + mq
. Accumulate mass
8:
Vp,sum ← Vp,sum + Vq
. Accumulate volume
9:mq ← NN+1 mq. Reduce mass
10:
Vq ← NN+1 Vq. Reduce volume
11:
end for12:
mp ← N1+1 mp,sum. Gathered mass
13:
Vp ← N1+1 Vp,sum. Gathered volume
14: end for
Algorithm 20 Compute New Field Values
1: for each node i in MPM grid do
2:vi ← 0. Initialize velocity
3:
Fi ← 0. Initialize deformation gradient4:
Jie ← 0. Initialize determinant5:bi ← 0. Initialize normalized strain
6:
wi,sum ← 0. Initialize weight sum
7: end for
8: for each particle p that existed before resampling do
9:
Mp ← all MPM grid nodes i satisfying wip > 0
10:for each MPM grid node i ∈ Mp do
11:w ← wip mp. Weight
12:
wi,sum ← wi,sum + w. Accumulate weight13:
v i ← v i + wv p
. Accumulate velocity
14:
Fi ← Fi + wFp . Accumulate deformation gradient15:
Ji ← Ji + wJp. Accumulate determinanteeebi ← bi + wbp. Accumulate normalized strain
16:
17:
end for
18: end for
19: for each node i in MPM grid do
20:
v i ← v i /wi,sum. Velocity
21:
Fi ← Fi /wi,sum. Deformation gradient22:
Ji ← Ji /wi,sum. Determinanteebi ← bi /wi,sum. Normalized strain
23:
24: end for
25: for each newly sampled particle p do
26:
vp ← 0. Initialize velocity
27:
Fp ← 0. Initialize deformation gradient28:
Jp ← 0. Initialize determinante29:
bp ← 0. Initialize normalized strain30:
Mp ← all MPM grid nodes i satisfying wip > 031:
for each MPM grid node i ∈ Mp do
32:v p ← v p + wip v i. Accumulate velocity
33:
Fp ← Fp + wip Fi. Accumulate deformation gradient34:
Jp ← Jp + wip Ji. Accumulate determinant
e
e
e
35:
bp ← bp + wip bi. Accumulate normalized strain
36:
end
for
e
e
e
37:
bp ← det(bp )−1/3 bp. Renormalize strain
1/3
38:
Fp ← Jp (det(Fp )−1/3 Fp ). Rescale deformation gradient39: end for

Algorithm 21 Merge Particles
1: P ← ∅. Set of potential particle pairs to merge
2: for each particle p do
3:
p.selected ← false
4: end for
5: for each particle p do
6:
Np ← all particles q satisfying |xq − xp | < 0.03r
7:
for each particle q do
8:
P ← P ∪ (p, q). Particles p and q are too close,. add the pair to the potential merge set
9:
end for
10: end for
11: D ← ∅. Set of particle pairs to merge
12: for each particle pair s ∈ P do
13:
if not s.p.selected and not s.q.selected then14:
D←s
. Append the particle pair15:s.p.selected ← true
. Skip s.p in future iterations
16:
s.q.selected ← true
. Skip s.q in future iterations
17:
end if
18: end for
19: for each particle pair s ∈ D do
20:m ← ms.p + ms.q
. Merge mass
21:
V ← Vs.p + Vs.q
. Merge volume22:v ← (ms.p v s.p + ms.q v s.q )/m. Merge velocity
23:
F ← (ms.p Fs.p + ms.q Fs.q )/m. Merge deformation gradient24:
J ← (ms.p Js.p + ms.q Js.q )/m. Merge determinanteeeb ← (ms.p bs.p + ms.q bs.q )/m. Merge strain
25:
e
e
e
26:
b ← det(b )−1/3 b. Renormalize strain
27:
F ← J 1/3 det(F)−1/3 F . Rescale deformation gradient28:
delete s.p and s.qe29:
add a particle with m, V , v, F, b
30: end for

ACM Transactions on Graphics, Vol. VV, No. N, Article XXX, Publication date: Month YYYY.