Efficient 2D Simulation on Moving 3D Surfaces
Dieter Morgenroth, Stefan Reinhardt, Daniel Weiskopf, Bernhard Eberhardt
AACM SIGGRAPH / Eurographics Symposium on Computer Animation 2020J. Bender and T. Popa(Guest Editors)
Volume 39 ( ), Number 8
Efficient 2D Simulation on Moving 3D Surfaces
D. Morgenroth , S.Reinhardt , , D. Weiskopf , and B. Eberhardt Media University Stuttgart, Germany Visualization Research Center, University of Stuttgart, Germany
Input Simulation Scalar Field Surface Simulation Rendering
Figure 1:
Water polluted with oil is poured into a cup. We simulate the oil film on top of the existing fluid simulation. The velocity and massare taken from the existing animation. On the left is the coarse input simulation. In the upper half of the middle image is the resulting scalarfield. The lower half shows the result of our method. Our high-resolution 2D simulation adds convincing visual details to the coarse inputsimulation. The right image displays the final rendering that needs the fine details from our surface simulation to generate the high-resolutionthin-film interference effects.
Abstract
We present a method to simulate fluid flow on evolving surfaces, e.g., an oil film on a water surface. Given an animatedsurface (e.g., extracted from a particle-based fluid simulation) in three-dimensional space, we add a second simulation on thisbase animation. In general, we solve a partial differential equation (PDE) on a level set surface obtained from the animatedinput surface. The properties of the input surface are transferred to a sparse volume data structure that is then used for thesimulation. We introduce one-way coupling strategies from input properties to our simulation and we add conservation of massand momentum to existing methods that solve a PDE in a narrow-band using the Closest Point Method. In this way, we efficientlycompute high-resolution 2D simulations on coarse input surfaces. Our approach helps visual effects creators easily integrate aworkflow to simulate material flow on evolving surfaces into their existing production pipeline.
CCS Concepts • Computing methodologies → Physical simulation;
1. Introduction
In typical workflows for generating digital visual effects, a team ofVFX artists iteratively refines a given sequence until they achievegood quality. Going from rough storyboards over blocked anima-tion to the final shot with many layers of physical simulations, eachshot passes through the VFX pipeline, where domain specialistsadd new effects and details. With physical simulations, often aneffect is finished and approved before secondary effects are lay-ered on top of it. For example, after simulating the fluid flow of awater surface, secondary effects like splashes and foam [TFK ∗ c (cid:13) (cid:13) a r X i v : . [ c s . G R ] S e p orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces Figure 2:
The seven steps of our method as a flow chart. The steps are placed in order of execution from left to right. After the input of theexternal data (Step 1), the data is transferred to the narrow-band (Step 2 and 3). Next, the surface is evolved (Step 4) and the outer processis coupled with the inner dynamics (Step 5), before solving the set of PDEs on the surface (Step 6). Finally, the fields are advected (Step 7). depicts the steps of our approach. We start with a moving surfaceas input geometry. After converting the input data into a distancefield and transferring values into a narrow-band grid around thesurface, we introduce quantities from the input 3D simulation intoour surface domain in a coupling step where the strength of thecoupling can be driven by parameters. Then, we use the ClosestPoint Method (CPM) to embed the 2-manifold in the 3D space ofthe moving surface and solve 2D PDEs in a 3D narrow-band.With our method, we can model secondary effects very effi-ciently, and we evaluate our approach by simulating a variety ofeffects such as pouring an oil film into water, simulating reaction-diffusion on a water surface, or the surfactants in a heated soapbubble. The source code is available on GitHub [MRWE20b].
2. Related Work
Secondary effects can be simulated on top of the base animation intwo different ways. The first way is to directly integrate them intothe simulation, for example, bubbles in foam [TSS ∗ ∗ ∗
15] on statictriangle meshes. Azencot et al. [AWO ∗
14] model fluids on trian-gle surfaces using their vorticity by a time-varying scalar function.Ruuth et al. [RM08] proposed the Closest Point Method (CPM)to calculate PDEs on surfaces. This was applied by Macdonald etal. [MMR13] to calculate reaction-diffusion terms on surfaces de-rived from point clouds. CPM on static surfaces can be calculatedin real time [AMT ∗ ∗
20] where the PDE issolved on a triangle mesh, or computing the flow within a soap bub-ble on a staggered spherical grid [HIK ∗ ∗
15] solved a PDE in a 3D narrow-band gridusing the CPM to enhance the original simulation with additionalturbulence. Closest to our approach is the Semi-Lagrangian Clos-est Point Method [AW13]. In comparison to existing CPM-basedmethods, we add equations for conservation of mass. This allowsus to correctly adapt scalar values when the surface area changes.We also adjust vector quantities to account for surface changes.Furthermore, we model a one-way coupling using mass and mo-mentum transfer to create a plausible linkage between the evolutionof the input surface and the behavior of the resulting 2D simulation.
3. Modeling PDEs on Evolving Surfaces
In this section, we describe the fundamentals of our method andhow we model the physics on the surface. We additionally describethe different physical phenomena that we model as secondary ef-fects on the surface. c (cid:13) (cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces Consider an evolving surface M . Such a surface can result from afinite-difference or Smoothed Particle Hydrodynamics (SPH) sim-ulation, a keyframe animation, or any time-dependent process. On M , we define another process (e.g., the simulation of a substance)with a set of PDEs and couple it to the base animation M resultsfrom. Therefore, the overall system is divided into three aspects: • Evolution of the surface M :Due to the surface evolution, the space on which we model thedynamic process is altered. In this step, we account for this fact,i.e., how quantities evolve alongside the surface. • Coupling :This aspect determines how the base animation affects the oneon the surface. • Modeling the dynamics on the surface :This step accounts for the secondary dynamic process and how itis modeled on the surface, i.e., how a set of PDEs can be solvedon the surface.An example of such a process could be a drop of ink on a plasticsheet moving in space. The first aspect (evolution of the surface)describes how the ink is transported in space when the plastic sheetis moving in normal direction. The second aspect (coupling) ac-counts for the movement of the plastic sheet in tangential direction,i.e., the acceleration of the ink due to friction forces, and the thirdaspect (the dynamic process on the surface) is concerned with thefluid behavior of the ink itself.
To model such a dynamic process on an evolving surface M ⊂ R ,we need to define scalar as well as vector-valued quantities on it.To this end, we attach scalar quantities a ( p , t ) ∈ R and vectorialquantities v : = v ( p , t ) ∈ T p M to every point p ∈ M , where T p M de-notes the tangent space at point p defined by the surface normal n : = n ( p , t ) and t denotes a certain point in time. Typical quantitieswould be, e.g., mass density and velocity for fluid flow on the sur-face. The velocity field u : = u ( p , t ) defines how the surface evolves,i.e., how M ( t + ∆ t ) results from M ( t ) , where ∆ t ∈ R > . For con-venience, we will write M instead of M ( t ) and M (cid:48) for M ( t + ∆ t ) .In the following, the surface evolution characterized by u will bereferred to as the outer process and the one on the surface will becalled inner dynamics.We construct a map O that describes the evolution of quantitiesdue to the outer process. To this end, we define the space M : M = (cid:91) p ∈ M { p } × R × T p M . (1)The quantities needed for the inner dynamics can now be describedas elements of M . Without loss of generality, we consider only onescalar and one vectorial quantity attached to point p . It is possibleto map an arbitrary number of quantities to p in the same way.The map O relates elements from M to elements in space M (cid:48) = (cid:83) p (cid:48) ∈ M (cid:48) { p (cid:48) } × R × T p (cid:48) M (cid:48) : O : M → M (cid:48) : p a v (cid:55)→ p (cid:48) a (cid:48) v (cid:48) . (2) Figure 3:
The map O relates p ∈ M and p (cid:48) ∈ M (cid:48) and also mapstangent space T p M into T p (cid:48) M (cid:48) . An illustration of O can be found in Fig. 3.We divide the influence of the outer process defined by u in nor-mal u n and tangential u t components. The tangential component isused to model friction between the outer and inner dynamics anddiscussed in Section 3.3. For now, we assume that the dynamics arecoupled without friction and, therefore, u t does not affect the innerdynamics, i.e., m ∈ M is only altered by u n .The rate of change D O Dt that a point p experiences from to actionof O is then defined by ˜ u n : = k u n , where D O Dt denotes the materialderivative and k ∈ R has to be chosen in such a way that p (cid:48) ∈ M (cid:48) holds. The motion of point p is directly governed by D O p Dt , in par-ticular its velocity is ˜ u n and, hence, D O p Dt = ˜ u n .A scalar quantity a is advected alongside the point p . We distin-guish between two types of scalar quantities: intensive and exten-sive ones [Red70]. An extensive property is a global property (e.g.,mass or volume). Such a property is only advected alongside p andnot changed, i.e., D O aDt =
0. In contrast, if a describes an intensivephysical property, such as density, we demand that0 = D O Dt (cid:90) M ( t ) a ( p , t ) dM ( p , t ) (3)holds true, where dM ( p , t ) are infinitesimal surface elements at po-sition p on the surface M at time t . In the following, we will skip t and p for convenience. Eq. 3 ensures that the total amount of sucha quantity does not change if no phenomena like mass transfer orchemical reactions are present. In this case, the rate of change of anintensive quantity a is given by D O aDt = − a ( ∇ · ( ˜ u n ) T p M ) , (4)where u T p M = ( I − nn T ) u is the velocity projected onto the tangentspace T p M at point p . This means that, if the surface diverges, theconcentration of a decreases, and vice versa. A detailed derivationis provided in Appendix A.When advecting a vectorial quantity v , we demand that v (cid:48) ∈ T p (cid:48) M (cid:48) holds after applying O . This can be ensured by considering v to be a small linear material line element subjected to the veloc-ity field ˜ u n . Its rate of change is the difference of the velocities atthe two ends of the element and can be described by D O v Dt = ∇ ˜ u n v .In our model, the vector field v describes the velocity of the innerdynamics. As the directions of the velocities are constrained by the c (cid:13) (cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces Figure 4:
The map O transforms velocity vector v to vector v (cid:48) .First, v is advected to the point p (cid:48) . Next, it is projected back intoT p (cid:48) M (cid:48) and scaled to conserve momentum. surface evolution, the total momentum will be altered and, there-fore, cannot be conserved. We can only ensure that the sum of theabsolute values does not change, i.e., we can enforce that0 = D O Dt (cid:90) M (cid:107) a v (cid:107) dM (5)will hold. Eq. 5 leads to D O (cid:107) v (cid:107) Dt = . (6)Hence, we need to adjust D O v Dt , so that the length of v is con-served. Using first-order Taylor expansion, the change of the lengthcan be computed by v · ∇ ˜ u n v v (cid:107) v (cid:107) . A detailed derivation can befound in Appendix B. This term is then subtracted to adhere Eq. 6.Metaphorically, this means that v is advected and rotated back ontothe surface. An illustration of the advection of a vectorial quantityis given in Fig. 4.If a is an intensive property, the combined material derivative ofthe map O then reads as D O Dt O = ˜ u n − a ( ∇ · ( ˜ u n ) T p M ) ∇ ˜ u n v − v · ∇ ˜ u n v v (cid:107) v (cid:107) . (7)Note that D O Dt a denotes the rate of change of a quantity a inducedby the operator O and D O Dt O denotes the combined material deriva-tive of O . So far we assumed that the outer process influences the inner dy-namics solely in normal direction. Next, we discuss how the tan-gential part u t = u − u n will influence the inner dynamics. Thisimplies that friction-less coupling is modeled, i.e., matter slides onthe surface and no adhesion is present. We model adhesive forcesto reduce the relative velocities v rel = v − u t between the outer pro-cess and inner dynamics as D c v Dt = − s v rel , (8)where s is a user-defined coefficient and D c Dt denotes the rate ofchange induced by coupling effects. Similar to Section 3.2, wecould also consider v to be a small linear material line elementexposed to the velocity u t and its rate of change would then read as (cid:0) ∇ T p M u t (cid:1) v , where ∇ T p M is defined as ∇ T p M = ( I − nn T ) ∇ . Wecombine both views and model the total rate of change of v causedby coupling effects as D c v Dt = − s v rel + s (cid:0) ∇ T p M u t (cid:1) v . (9)The outer process induces sinks and sources to the inner dynam-ics modeled on the surface. Therefore, we model transfer of in-tensive physical quantities (e.g., density) from the outer process toinner dynamics via sinks and sources on the surface. Inspired byFick’s laws of diffusion, we want that the concentration a out fromthe outer process and the concentration a of the inner dynamicswill align. In other words, the concentration of a substance on thesurface will be changed by D c aDt = − s a rel , (10)to reduce the concentration difference a rel = a − a out . The speedof the reduction is determined by the factor s . If only sources aremodeled, the concentration difference is restricted to negative val-ues, i.e., a rel ≤
0. Sinks are modeled by restricting a rel ≥ After discussing the surface evolution and coupling, we describehow dynamics on the surface can be modeled. To describe such aprocess we use a set of PDEs, e.g., D s v Dt = F ( t , a , v , ∂ a , ∂ a , ..., ∂ v , ∂ v , ... ) and (11) D s aDt = F ( t , a , v , ∂ a , ∂ a , ..., ∂ v , ∂ v , ... ) , (12)where the functions F and F characterize the phenomena mod-eled on the surface. The term ∂ i : = (cid:26) ∂ | α | ∂ x α ∂ x α ∂ x α : | α | = i (cid:27) de-notes the set of all partial derivatives of the order i , with α beinga multi-index. For example, if fluid flow is modeled, F describesthe momentum conservation and F the mass conservation of theNavier-Stokes equations.To solve Eqs. 11 and 12 on the surface, we only allow for tan-gential deviations. To this end, we solve them locally in the corre-sponding tangent spaces T p M . This implies that, instead of usingthe ordinary spatial derivations, we project them onto the surface,e.g., the operator ∇ is replaced by ∇ T p M . The governing equationsfor our model are then given by: D v Dt = D O v Dt + D c v Dt + D s v Dt and (13) D aDt = D O aDt + D c aDt + D s aDt . (14)While the outer process is ignored in Equations 11 and 12, it isincluded in Equations 13 and 14, i.e., they govern the overall dy-namics on the surface.Next, we describe different physical phenomena that we use inthis paper to show the versatility of our model. c (cid:13) (cid:13)(cid:13)
0. Sinks are modeled by restricting a rel ≥ After discussing the surface evolution and coupling, we describehow dynamics on the surface can be modeled. To describe such aprocess we use a set of PDEs, e.g., D s v Dt = F ( t , a , v , ∂ a , ∂ a , ..., ∂ v , ∂ v , ... ) and (11) D s aDt = F ( t , a , v , ∂ a , ∂ a , ..., ∂ v , ∂ v , ... ) , (12)where the functions F and F characterize the phenomena mod-eled on the surface. The term ∂ i : = (cid:26) ∂ | α | ∂ x α ∂ x α ∂ x α : | α | = i (cid:27) de-notes the set of all partial derivatives of the order i , with α beinga multi-index. For example, if fluid flow is modeled, F describesthe momentum conservation and F the mass conservation of theNavier-Stokes equations.To solve Eqs. 11 and 12 on the surface, we only allow for tan-gential deviations. To this end, we solve them locally in the corre-sponding tangent spaces T p M . This implies that, instead of usingthe ordinary spatial derivations, we project them onto the surface,e.g., the operator ∇ is replaced by ∇ T p M . The governing equationsfor our model are then given by: D v Dt = D O v Dt + D c v Dt + D s v Dt and (13) D aDt = D O aDt + D c aDt + D s aDt . (14)While the outer process is ignored in Equations 11 and 12, it isincluded in Equations 13 and 14, i.e., they govern the overall dy-namics on the surface.Next, we describe different physical phenomena that we use inthis paper to show the versatility of our model. c (cid:13) (cid:13)(cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces Figure 5:
A rotating sphere where we couple the velocity of the in-ternal simulation with the velocity of the sphere. Details are addedwith artificial vorticity and coupling noise. Simulation resolution is207 × × We model different physical phenomena on the surface and, there-fore, employ different sets of equations.
Fluid Flow.
One application of our model is the simulation of fluidflow on an evolving surface. Figs. 5, 7, and 8b show examples ofsuch a flow. The surface rotates and, due to friction forces, the fluidon the surface starts to move. To model viscous fluid flow on thesurface we use the Navier-Stokes momentum equation: D s v Dt = − ρ ∇ T p M P + ν ∇ T p M v + ρ F b , (15)where v is the fluid’s velocity, ν the viscosity constant, P is thepressure, and F b are body forces (e.g., gravity).The evolution of density ρ is modeled by using the general con-tinuity equation: D s ρ Dt = σ − ρ ( ∇ T p M · v ) , (16)where σ is the rate of generation of the substance per unit volume,i.e., it is used to model sinks ( σ <
0) or sources ( σ > F ( ρ , v ) =
0. As aresult, Eq. 16 simplifies to a volume conservation law: ∇ T p M · v = σρ . (17)If no sinks or sources are present (i.e., s = ∇ T p M · v =
0. Note that if the outer process is divergence-free, i.e., ∇ · u =
0, we obtain an incompressible flow because D ρ Dt = Buoyancy-induced Flow.
Fig. 6 shows a static hemisphere. Someareas on the surface are heated, others cooled by outer tempera-ture constraints. We model temperature transfer according to Eq. 10and, therefore, the fluid on the surface changes its temperature. Nat-ural convection arises due to temperature differences in the fluid.Such buoyancy-driven flow can be modeled using the so-called
Figure 6:
Thermal convection on a hemisphere. The temperatureis color-coded, where the temperature rises from blue over whiteto red. Some areas are heated and others are cooled. Due to thesetemperature differences, buoyancy-driven flow arises. Simulationresolution is 207 × × Boussinesq approximation [Bou97]. The Boussinesq approxima-tion assumes incompressible flow and that variations of densityonly occur due to temperature differences: ρ = ρ − βρ ( T − T ) ,where β is the coefficient of thermal expansion, T the referencetemperature, and ρ the reference density. We assume an ideal gasand, therefore, β = T . If we assume gravity as the only body-forcepresent, Eq. 15 becomes D s v Dt = − ρ ∇ T p M P + ν ∇ T p M v + (cid:18) − TT (cid:19) g , (18)with g being the gravitational constant. To model changes in thetemperature field we use the convection-diffusion equation: D s T Dt = µ d ∇ T p M T , (19)where we assume a constant diffusion coefficient µ d . Reaction-diffusion.
Reaction-diffusion is commonly used tomodel chemical reactions of one or more substances (e.g., a sub-stance is transformed into another due to chemical reactions) andthe diffusion of the substance(s) in space. Fig. 8a shows such a pro-cess. Two chemical substances diffuse and react with each other. Inthis example, we model a two-component reaction-diffusion pro-cess given by D s f Dt = R ( f , f ) + µ d ∇ T p M f and (20) D s f Dt = R ( f , f ) + µ d ∇ T p M f . (21)The functions R and R are the reaction terms and characterize thesystem.We use the Gray Scott Model [GS84] to define the reactionterms, i.e., to model how different morphogens react: R ( f , f ) = − f f + β ( − f ) , (22) R ( f , f ) = f f − ( β + γ ) f , (23)To generate the example illustrated in Fig. 8a, we set the parametersto β = . γ = . µ d = .
1, and µ d = . c (cid:13) (cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces Figure 7:
Fluid flow simulation on top of a FLIP fluid riverbed sim-ulation. Performing a surface simulation of a color-band demon-strates that we can add fine-scaled details. The split screen showsthe input fluid simulation on the left and the added surface flow onthe right.
4. Connecting to the Base Animation
In this section, we describe how we apply the model from Section 3to create a simulation system. Our method can be divided into sevensteps, as illustrated in Fig. 2. We start with a coarse input (Fig. 2,Step 1), e.g., a fluid simulation. Then, we create a signed distancefield from this simulation (Fig. 2, Step 2) and write the simulationproperties such as velocity, density, or color into 3D fields that arelaid out in a narrow-band grid, which is our generic input format(Fig. 2, Step 3). We bring quantities from the outer process intoour surface domain in a coupling step and compensate for the sur-face evolution as described in Section 4.4 (Fig. 2, Step 4 and 5).Then, we simulate 2D details to enhance the coarse input by solv-ing a PDE in the 2D surface space (Fig. 2, Step 6) that results in avelocity field that advects all grids in the last step (Fig. 2, Step 7).To solve a PDE on a surface we use the Closest Point Method(CPM) [RM08] calculated in a 3D narrow-band. The main idea ofthe CPM is that if the values of a field are constant along the normalof the surface, many equations calculated in this 3D field are thesame as if they were calculated in the tangent spaces of the surface.Section 4.3 describes the process of converting 3D fields into fieldsthat have constant values along the surface normal.We implemented our method where each part is interchangeableusing a system of modules that change data flowing through thesystem. In the following, we will describe each step in detail.
The first step implements data input and is able to process a varietyof base animation types. In our examples, we consider keyframeanimation (Fig. 5), SPH particle simulation (Fig. 12 and 8a), andfluid implicit particle (FLIP) simulation (Fig. 7). The only require-ment for an input is that we can convert the surface geometry datainto a signed distance field and that values for the surface velocitycan be generated on the surface. (a)
Reaction-diffusion (b)
Fluid flow
Figure 8:
With our approach, we can model different behaviors onthe same input data. The upper image shows a reaction-diffusionsimulation on top of a dam break SPH simulation. The lower imageshows a fluid flow on top of the same input simulation. Simulationresolution is 199 × × For converting polygon meshes to signed distance fields, severalmethods are available, e.g., [SGGM06], [XB14], or [KDBB17]. Aspecial case arises when we convert particle systems to distancefields. Inside the fluid, there are particles everywhere. Here, thedistance field obtained by these methods is not usable since weneed the distance to the surface and not to the nearest particle.For particle-based simulations, we can either transform the sim-ulation to a polygon mesh first or go directly from particles tosigned distance fields by defining an implicit surface from the par-ticles [ZB05] and then take the distance to this implicit surface.
We create a narrow-band around the surfaces to capture the veloc-ity, density, and other properties of the flow from the nearest surfacepoint. We only fill cells near the surface based on the distance field.Cells that a farther away than a certain threshold are left empty,as illustrated in Fig. 9. We keep track of two velocity fields, theinternal velocity field (denoted as v in Section 3) and the externalvelocity field (denoted as u ). We write the velocity of the incomingsurface into the external velocity grid. Depending on the use case,we either initialize the internal velocity grid with a snapshot of theexternal velocity (for example, Fig. 12), or we create a custom ini-tialization routine to define the initial internal velocity, for example,by initializing with vanishing velocity (Fig. 5). We also keep track c (cid:13) (cid:13)(cid:13)
We create a narrow-band around the surfaces to capture the veloc-ity, density, and other properties of the flow from the nearest surfacepoint. We only fill cells near the surface based on the distance field.Cells that a farther away than a certain threshold are left empty,as illustrated in Fig. 9. We keep track of two velocity fields, theinternal velocity field (denoted as v in Section 3) and the externalvelocity field (denoted as u ). We write the velocity of the incomingsurface into the external velocity grid. Depending on the use case,we either initialize the internal velocity grid with a snapshot of theexternal velocity (for example, Fig. 12), or we create a custom ini-tialization routine to define the initial internal velocity, for example,by initializing with vanishing velocity (Fig. 5). We also keep track c (cid:13) (cid:13)(cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces Figure 9:
The distance field is stored in a narrow-band around thesurface. In the same way, velocity and scalar fields are stored in asparse data structure. of the scalar properties that we use in our simulation, e.g., the con-centration of substances for the reaction-diffusion example.To use the CPM our inputs must be narrow-band grids wherethe values are constant along the normal direction of the surface.Using an integer look-up grid that stores the cell coordinates of theclosest point of the surface in each grid cell, we can fill in a gridwith values from the closest surface point. We refer to this step asthe CPM extension. Here, we extend the values near the surfacealong the normal direction. The CPM extension is our final step inthe narrow-band creation phase.
Before we use the velocities or scalar values like color or density ina PDE, we first execute the correction steps for density and velocity,as described in Section 3.For mass conservation, we have to solve Eq. 4. As a buildingblock for this calculation, we need a module that implements theoperator ∇ · ( ˜ u n ) T p M . This operator will calculate a scalar value forthe divergence of the velocity field but with respect to the tangentspace of a surface defined by the gradient of the input distance field.The result is then applied to the scalar fields that need correction.The interesting part of this divergence calculation is that we usethe original 3D velocities for the divergence operator. We need totake the velocities of adjacent cells into account but projected tothe surface normal of the current grid value, not the normals of therespective neighbor cells. An example to better understand the dif-ference is the one of a growing sphere. Although the 2D velocitiesin surface space are zero at each point, the divergence is not. To getcorrect 2D divergence where a growing sphere causes sinks, and ashrinking sphere creates sources in the mass conservation equation,we have to project the velocities of neighbor cells using the normalof the current cell.For the conservation of momentum, the “Project Vector” modulealters the velocities based on the surface tangent. Given a sourcevector grid and a distance field gradient, this module will projectthe input vectors onto the surface by subtracting the normal vectorcomponent from the gradient field. There is an option to maintainthe length from the input vector, which resembles rotating the vec-tor down onto the surface, as shown in Fig. 4. To apply Eq. 9, we Time M ea n D e n s it y SimulationAnalytic
Figure 10:
Simulated density on a growing and shrinking spherecompared to the analytic solution. With our simulation, we can re-produce the analytical result. The numerical diffusion is negligible. add a module to calculate the jacobian of a velocity field and useit to change the velocities of a second velocity field. Both modulesare applied to the velocity field for the internal velocity.After the values are adapted based on the surface evolution, wecouple them to the values from the underlying simulation.
For coupling the velocities, we send the velocity grids to a modulethat applies Eq. 9 to the inner velocity. The parameters for s s s As mentioned, we are operating on 3D grids where values do notchange along the normal direction. The CPM allows us to solvePDEs in the surface space using these grids. We implemented dif-ferent sets of PDEs that can all be solved using our new modules.The modules work in 3D, but we specifically designed them forthe fields that result from the CPM extension. Due to the CPM, thegradient naturally operates in surface direction, but the divergence,curl, and Jacobian have to be changed: When applying them, weuse the projected version as described in Section 3, i.e., the inputvector is projected onto the surface.To simulate fluid flow on the surface, we solve Eq. 15.Here, we use the projected divergence operator ∇ · ( ˜ u n ) T p M . The“Divergence-Free” module will remove divergence with respect tothe tangent space of a surface and take the divergence of the surfaceevolution into account. For the buoyancy-induced flow, we scaledthe gravity depending on temperature as described in Eq. 18 andadded a diffusion module for temperature diffusion. To simulate c (cid:13) (cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces reaction-diffusion for the example shown in Fig. 8a we wrote amodule to implement Equations 20–23. As a last step, we advectall fields using Eq. 13, i.e., with the resulting velocities. Our modular approach can be easily implemented with existingframeworks like the SideFX TM Houdini software package, whichis widely used for VFX creation and offers an integration ofthe OpenVDB libraries [MLJ ∗
13] as a sparse volume represen-tation for the calculations. We separated our algorithm into in-dependent modules and implemented them as a separate opera-tors, to be used in Houdini’s simulation node graph framework.The OpenVDB framework offers a data structure, the OpenVDBgrid, to create narrow-bands. All subsequent calculations employthis structure and only spend computation time where needed. Our“CPM Extension” module is the central block to build CPM so-lutions in Houdini. Here, we implemented a nearest-neighbor in-terpolation as suggested by Kim et al. [KTT13], but also box andquadratic interpolation. To generate the SPH base simulations, weused divergence-free SPH [BK15] with consistent Shepard interpo-lation [RKEW19]. For more details about SPH, we refer the readerto the SPH tutorial by Koschier et al. [KBST19]. When simulatingfluid flow, we used vorticity confinement [FSJ01].
5. Results5.1. Versatility and Simulation Quality
To show the versatility of our method, we modeled different phys-ical phenomena as described in Section 3.5, and tested them on avariety of scenarios. We chose base animations with different char-acteristics, from static (Fig. 6) to hand-animated meshes (Fig. 5),from coarse-scale SPH-based fluid simulations, including suddenchanges in topology (Fig. 8), to high-resolution multi-phase SPHsimulations (Fig. 12). In all of the above situations, we were ableto add a fine-scale secondary simulation on top of these surfaces,revealing fine-scale details as promised.With our approach, we can add effects onto the surface of thesimulation. Simulating a second fluid flow on the surface enablesone to add another phenomena on top of an existing simulation. Asshown in the oil film example (Fig. 12), the effect of oil spreadingon the surface is achieved by simulating a second, fine-scale fluidflow on the surface, which is just added to the base-simulation. Asmentioned, the level of detail of the secondary simulation can bechosen independently and increasing the simulation resolution ofthe underlying SPH simulation would not have the same effect. In-stead, we would just have a scalar field with higher resolution, likethe one shown in Fig. 12b, but would not have simulated the lawsof the secondary flow. Also, due to the modeled mass transfer, theoil particles that emerge from below the surface can contribute toour 2D simulation. This coupling introduces significantly more de-tails than just a plain simulation on the surface itself. The amount ofdetail added is significant even for low-resolution base animations,and it can be further improved (Fig. 12c to 12d) by increasing theresolution. Fig. 8a presents a reaction-diffusion of two chemicalssimulated on top of a coarse dam break simulation. This exampleillustrates that we can simulate not only a second flow equation on .
01 0 .
019 0 .
028 0 .
037 0 .
046 0 .
055 0 .
065 0 .
074 0 . Gridsize M ea n A b s o l u t e E rr o r Figure 11:
The mean error for different grid sizes. The smaller thegrid cells, the closer the values are to the analytic solution. the surface but any kind of desired physical phenomena that can bedescribed by a set of PDEs.In the buoyancy-induced flow example (Fig. 6), we simulatethermal convection and demonstrate emergence of isolated vorticeson a static surface, replicating the physical experimental setup bySeychelles et al. [SABK08]. The setup consists of a static hemi-sphere on a heating plate giving rise to thermal convection. Wewere able to create fine detailed vortex structures using the Boussi-nesq approximation.A hand-animated sphere with a periodically oscillating radiuswas used to test our approach on mass conservation. The meandensity per surfactant was calculated and plotted as a function overtime as a graph (Fig. 10). In this example, it is possible to analyt-ically calculate the density change that is needed to ensure massconservation and compare it with our simulated result. As shownin Fig. 10, we can reproduce the analytical solution over time. Thesmall loss in density can be attributed to numerical diffusion and isnegligible in this case. To test the accuracy depending on the gridresolution, we simulated this scenario with different grid resolu-tions. As can be seen in Fig. 11, our mass conservation approachworks as expected and the simulation converges to the analyticalsolution using finer grid resolutions. Further investigations on theconvergence order is left for future work.We refer the reader to the accompanying video to watch some ofthe examples in motion. Supplementary material discusses the vec-tor length conservation experiment. Moreover, we provide a pre-compiled version of our plugin for Houdini with two sample scenesfor testing on Zenodo [MRWE20a].
The performance was measured for every operator individuallywith the rotating sphere example (Fig. 5). The grid resolution hasthe most significant impact on the computation time, but also thewidth of the narrow-band profoundly influences the performance.For the rotating sphere example, the overall computation time perframe with 5 substeps was 28 s. The soap film example from Fig. 12took approximately 28 s per frame with 2 substeps. The measure-ments were taken on a notebook with an Intel(R) Core i7-6700KCPU and a NVidia GeForce GTX 1070. The simulation of theunderlying fluid is not included in the measurements. As can be c (cid:13) (cid:13)(cid:13)
The performance was measured for every operator individuallywith the rotating sphere example (Fig. 5). The grid resolution hasthe most significant impact on the computation time, but also thewidth of the narrow-band profoundly influences the performance.For the rotating sphere example, the overall computation time perframe with 5 substeps was 28 s. The soap film example from Fig. 12took approximately 28 s per frame with 2 substeps. The measure-ments were taken on a notebook with an Intel(R) Core i7-6700KCPU and a NVidia GeForce GTX 1070. The simulation of theunderlying fluid is not included in the measurements. As can be c (cid:13) (cid:13)(cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces (a) Source SPH particles (b)
Original scalar field derived fromparticle attributes (c)
2D simulation with grid resolutionof 0.05 (d)
2D simulation with grid resolutionof 0.02
Figure 12:
Pouring polluted water into a bowl. The fine-grained simulation enhances the underlying SPH simulation. Our mass transfereven captures oil particles that emerge from under the surface. Different resolutions can be generated independent from the input resolution.This allows iterative refinement of parameters as needed in typical VFX production workflows.
Narrow BandPDECouplingSurface EvolutionAdvectionOther Divergence FreeJacobianCPM Extension ICPM Extension IICPM Extension IIICPM Extension IVDistance FieldVDB OperationsVorticityAdhesive ForcesVDB OperationsProject Vector IProject Vector IIDisplay Compute Outer VelocityInitial Mesh Generation
Figure 13:
The calculation costs of our nodes for the CPM meth-ods are of the same magnitude as the existing Houdini operatorstypically used in simulations. seen in Fig. 13, most of the computation time for a single task isspent in the “Divergence Free” operator. This node uses the Open-VDBs Poisson solver with a maximum of 100 iterations to obtaina divergence-free vector field. Second most computation time isspent to compute the Jacobian for the coupling step. In third placeis Houdini’s advection operator. All other operators have smallercomputation costs. This shows that the computation time for oursteps is of the magnitude as Houdini’s advection node.
6. Conclusion and Future Work
We presented a method to solve PDEs on evolving surfaces wherewe considered the external movement as the driving process.We derived physically motivated conservation laws and couplingstrategies. We coupled simulations of 3D and 2D space in a waythat allows us to compute high-resolution 2D simulations on coarseinput surfaces efficiently and that exposes physically plausible pa-rameters to control the strength of the coupling. We demonstratedthat our approach can be used for different types of problems thatrequire solving PDEs on a surface. We showed how to integrate all necessary steps into a VFX production environment in a modularway so that the building blocks can be reused and we provide areference implementation as open source. The 2D surface simula-tion showed to be a usable tool for VFX creation that can add abig visual difference to scenes. By using the CPM there are somelimitations to our approach. As most CPM implementations, ourmethod cannot model surfaces with hard edges or high frequencies.These limitations restrict use-cases to smooth surfaces. So far, weassume a surface without boundaries. Future research could inte-grate different boundary conditions. Testing further types of PDEsto achieve effects like droplet formation or fingering as presentedby Vantzos et al. [AVW ∗
15] but on moving surfaces is another topicfor future work. Another area of interest is reinjecting 3D fluidsfrom the surface as a simple form of two-way coupling of the sim-ulations, e.g., for water droplets.
Acknowledgements
We want to thank Mariusz Wesierski for his help in setting up theriver simulation and his support in all Houdini questions. This workis partly supported by “Kooperatives Promotionskolleg Digital Me-dia” at Stuttgart Media University and the University of Stuttgart.
Appendix A: Conservation of Intensive Scalar Quantities
In Section 3.2, we stated that Eq. 4 directly follows from Eq. 3if a represents an intensive scalar quantity. In the following, weprovide a detailed derivation of this fact. We want to ensure that thetotal amount of a does not change if there are no phenomena likemass transfer or chemical reactions present. As the surface evolves,the surface density a needs to be adjusted to ensure that the totalamount stays constant. Changing the differentiation and integrationorder in Eq. 3 results in0 = (cid:90) M D O Dt ( a dM ) = (cid:90) M (cid:18) D O aDt dM + a D O dMDt (cid:19) . (24)Since dM evolves, we need to evaluate the term D O dMDt .According to Stone [Sto90], the material derivative of a surface c (cid:13) (cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces element can be described by DDt dM = ( ∇ T p M · ˜ u n ) dM , (25)where ∇ T p M = ( I − nn T ) ∇ is the gradient operator projected ontothe tangential space T p M at point p , and ˜ u n is the velocity evolvingthe surface. Without loss of generality, we assume that the outerprocess is described by u instead of ˜ u n and use DdMDt instead of D O dMDt to derive the rate of change of dM . To describe Eq. 25 inmore detail, we first show that the identity DDt d M = ( ∇ · u ) d M − ( ∇ u ) T d M (26)holds for an arbitrary oriented surface element d M = n dM [Bat00].To this end, we first investigate how volume elements dV and lineelements d l are changed due to the base animation. We consider theelements to be that small, that they are subject only to pure strainingmotion and rigid rotations [Bat00].The rate of change of the volume element dV is described by D dVDt = ∇ · u dV . (27)We assume that the line element d l is linear (and, therefore, canbe described by a vector) and stays approximately straight. Underthese assumptions, its rate of change is simply the difference of thevelocities at the two ends of the element and can be described by D d l Dt = ∇ u d l . (28)The volume element dV can also be described by dV = d M · d l .Inserting dV into Eq. 27 results in: ( ∇ · u ) d M · d l = ∇ · u dV = D ( d M · d l ) Dt = d M · D d l Dt + D d M Dt · d l = d M · ( ∇ u d l ) + D d M Dt · d l = (( ∇ u ) T d M ) · d l + D d M Dt · d l . (29)As Eq. 29 has to hold for arbitrary line elements d l , we obtainEq. 26. To get the formulation in Eq. 25 we take the inner-productof n with Eq. 26 and use the identity n T ( ∇ u ) T n = ( nn T ∇ ) · u : DDt dM = DDt n T n dM = n T DDt d M = n T ( ∇ · u ) d M − n T ( ∇ u ) T d M = ( ∇ · u ) dM − ( nn T ∇ ) · u dM = (( I − nn T ) ∇ ) · u dM = ( ∇ T p M · u ) dM . (30)Inserting Eq. 25 into Eq. 24 results in0 = (cid:90) M ( t ) (cid:18) D O aDt dM + a ( ∇ T p M · ˜ u n ) dM (cid:19) = (cid:90) M ( t ) (cid:18) D O aDt + a ( ∇ T p M · ˜ u n ) (cid:19) dM . (31) As Eq. 31 holds for an arbitrary surface M , we get0 = D O aDt + a ( ∇ T p M · ˜ u n ) = D O aDt + a ( ∇ · ( ˜ u n ) T p M ) . (32)Therefore, the rate of change of a intensive quantity a is defined by D O aDt = − a ( ∇ · ( ˜ u n ) T p M ) . (33) Appendix B: Length-preserving Evolution of a Vector Field
In Section 3.2, we stated that Eq. 5 is satisfied when setting D O v Dt = ∇ ˜ u n v − v · ∇ ˜ u n v v (cid:107) v (cid:107) . (34)To confirm that Eq. 34 results from Eq. 5, we first show that D O (cid:107) v (cid:107) Dt = = D O Dt (cid:90) M (cid:107) a v (cid:107) dM = D O Dt (cid:90) M a (cid:107) v (cid:107) dM = (cid:90) M (cid:18) D O aDt (cid:107) v (cid:107) dM + a D O (cid:107) v (cid:107) Dt dM + a (cid:107) v (cid:107) D O dMDt (cid:19) = (cid:90) M a D O (cid:107) v (cid:107) Dt dM . (35)This is true for an arbitrary surface M and we obtain: D O (cid:107) v (cid:107) Dt = . (36)To achieve this, we evolve v with ∇ ˜ u n and compensate lengthchanges. The change in length of v can be expressed by v · ∇ ˜ u n v (cid:107) v (cid:107) v (cid:107) v (cid:107) , (37)and we obtain: D O v Dt = ∇ ˜ u n v − ( v · ∇ ˜ u n v ) v (cid:107) v (cid:107) . (38)This can be shown when writing the derivative as the limit of thedifference quotients. To this end, we use the first-order Taylor ex-pansion of v with respect to t to define ˆ v = v + ∆ t ∇ ˜ u n v . Compen-sating for length changes, v (cid:48) : = v ( t + ∆ t ) can be approximated by: v (cid:48) = ˆ v − ( (cid:107) ˆ v (cid:107) − (cid:107) v (cid:107) ) ˆ v (cid:107) ˆ v (cid:107) (39) c (cid:13) (cid:13)(cid:13)
In Section 3.2, we stated that Eq. 5 is satisfied when setting D O v Dt = ∇ ˜ u n v − v · ∇ ˜ u n v v (cid:107) v (cid:107) . (34)To confirm that Eq. 34 results from Eq. 5, we first show that D O (cid:107) v (cid:107) Dt = = D O Dt (cid:90) M (cid:107) a v (cid:107) dM = D O Dt (cid:90) M a (cid:107) v (cid:107) dM = (cid:90) M (cid:18) D O aDt (cid:107) v (cid:107) dM + a D O (cid:107) v (cid:107) Dt dM + a (cid:107) v (cid:107) D O dMDt (cid:19) = (cid:90) M a D O (cid:107) v (cid:107) Dt dM . (35)This is true for an arbitrary surface M and we obtain: D O (cid:107) v (cid:107) Dt = . (36)To achieve this, we evolve v with ∇ ˜ u n and compensate lengthchanges. The change in length of v can be expressed by v · ∇ ˜ u n v (cid:107) v (cid:107) v (cid:107) v (cid:107) , (37)and we obtain: D O v Dt = ∇ ˜ u n v − ( v · ∇ ˜ u n v ) v (cid:107) v (cid:107) . (38)This can be shown when writing the derivative as the limit of thedifference quotients. To this end, we use the first-order Taylor ex-pansion of v with respect to t to define ˆ v = v + ∆ t ∇ ˜ u n v . Compen-sating for length changes, v (cid:48) : = v ( t + ∆ t ) can be approximated by: v (cid:48) = ˆ v − ( (cid:107) ˆ v (cid:107) − (cid:107) v (cid:107) ) ˆ v (cid:107) ˆ v (cid:107) (39) c (cid:13) (cid:13)(cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces and the derivative of v is written as: D O v Dt = lim ∆ t → v (cid:48) − v ∆ t = lim ∆ t → ∆ t (cid:18) ∆ t ∇ ˜ u n v − ( (cid:107) ˆ v (cid:107) − (cid:107) v (cid:107) ) ˆ v (cid:107) ˆ v (cid:107) (cid:19) = lim ∆ t → (cid:18) ∇ ˜ u n v − ∆ t (cid:18) − (cid:107) v (cid:107)(cid:107) ˆ v (cid:107) (cid:19) v + (cid:18) − (cid:107) v (cid:107)(cid:107) ˆ v (cid:107) (cid:19) ∇ ˜ u n v (cid:19) = lim ∆ t → (cid:32) ∇ ˜ u n v − (cid:107) ˆ v (cid:107) − (cid:107) v (cid:107) ∆ t (cid:107) ˆ v (cid:107) ( (cid:107) ˆ v (cid:107) + (cid:107) v (cid:107) ) v + (cid:18) − (cid:107) v (cid:107)(cid:107) ˆ v (cid:107) (cid:19) ∇ ˜ u n v (cid:19) = lim ∆ t → (cid:32) ∇ ˜ u n v − ∆ t v · ( ∇ ˜ u n v ) + ∆ t (cid:107)∇ ˜ u n v (cid:107) ∆ t (cid:107) ˆ v (cid:107) ( (cid:107) ˆ v (cid:107) + (cid:107) v (cid:107) ) v + (cid:18) − (cid:107) v (cid:107)(cid:107) ˆ v (cid:107) (cid:19) ∇ ˜ u n v (cid:19) = lim ∆ t → (cid:32) ∇ ˜ u n v − v · ( ∇ ˜ u n v ) + ∆ t (cid:107)∇ ˜ u n v (cid:107) (cid:107) ˆ v (cid:107) ( (cid:107) ˆ v (cid:107) + (cid:107) v (cid:107) ) v + (cid:18) − (cid:107) v (cid:107)(cid:107) ˆ v (cid:107) (cid:19) ∇ ˜ u n v (cid:19) = ∇ ˜ u n v − v · ( ∇ ˜ u n v ) (cid:107) v (cid:107) v (40) References [ADAT13] A
KINCI
N., D
IPPEL
A., A
KINCI
G., T
ESCHNER
M.: Screenspace foam rendering.
Journal of WSCG 21 , 3 (2013), 173–182. 2[AMT ∗
12] A
UER
S., M
ACDONALD
C., T
REIB
M., S
CHNEIDER
J.,W
ESTERMANN
R.: Real-time fluid effects on surfaces using the clos-est point method.
Computer Graphics Forum 31 , 6 (2012), 1909–1923. doi:10.1111/j.1467-8659.2012.03071.x . 2[AVW ∗
15] A
ZENCOT
O., V
ANTZOS
O., W
ARDETZKY
M., R
UMPF
M.,B EN -C HEN
M.: Functional thin films on surfaces. In
Proceedings of the14th ACM SIGGRAPH / Eurographics Symposium on Computer Anima-tion (2015), pp. 137–146. doi:10.1145/2786784.2786793 . 2,9[AW13] A
UER
S., W
ESTERMANN
R.: A semi-Lagrangian closest pointmethod for deforming surfaces.
Computer Graphics Forum 32 , 7 (2013),207–214. doi:10.1111/cgf.12228 . 2[AWO ∗
14] A
ZENCOT
O., W
EISSMANN
S., O
VSJANIKOV
M.,W
ARDETZKY
M., B EN -C HEN
M.: Functional fluids on sur-faces.
Computer Graphics Forum 33 , 5 (2014), 237–246. doi:10.1111/cgf.12449 . 2[Bat00] B
ATCHELOR
G. K.:
An Introduction to Fluid Dynamics , 8th ed.Cambridge Mathematical Library. Cambridge University Press, 2000. doi:10.1017/CBO9780511800955 . 10[BK15] B
ENDER
J., K
OSCHIER
D.: Divergence-free smoothed par-ticle hydrodynamics. In
Proceedings of the 14th ACM SIG-GRAPH/Eurographics Symposium on Computer Animation (2015),pp. 147–155. doi:10.1145/2786784.2786796 . 8[Bou97] B
OUSSINESQ
J.:
Théorie de l’écoulement tourbillonnant et tu-multueux des liquides dans les lits rectilignes a grande section , 1 ed.Gauthier-Villars, 1897. 5[CPPK07] C
LEARY
P. W., P YO S. H., P
RAKASH
M., K OO B. K.: Bub-bling and frothing liquids.
ACM Transactions on Graphics 26 , 3 (2007),97:2–97:6. doi:10.1145/1276377.1276499 . 2[FSJ01] F
EDKIW
R., S
TAM
J., J
ENSEN
H. W.: Visual simulation ofsmoke. In
Proceedings of the 28th Annual Conference on ComputerGraphics and Interactive Techniques (2001), pp. 15–22. doi:10.1145/383259.383260 . 8[GDP16] G
AGNON
J., D
AGENAIS
F., P
AQUETTE
E.: Dynamic lappedtexture for fluid simulations.
The Visual Computer 32 , 6 (2016), 901–909. doi:10.1007/s00371-016-1262-8 . 2 [GS84] G
RAY
P., S
COTT
S.: Autocatalytic reactions in the isothermal,continuous stirred tank reactor: Oscillations and instabilities in the sys-tem a + 2b → → c. Chemical Engineering Science 39 , 6 (1984),1087–1097. doi:10.1016/0009-2509(84)87017-7 . 5[HH16] H
ILL
D. J., H
ENDERSON
R. D.: Efficient fluid simulation onthe surface of a sphere.
ACM Transactions on Graphics 35 , 2 (2016),16:1–16:9. doi:10.1145/2879177 . 2[HIK ∗
20] H
UANG
W., I
SERINGHAUSEN
J., K
NEIPHOF
T., Q U Z.,J
IANG
C., H
ULLIN
M. B.: Chemomechanical simulation of soap filmflow on spherical bubbles.
ACM Transactions on Graphics 39 , 4 (2020). doi:10.1145/3386569.3392094 . 2[IAAT12] I
HMSEN
M., A
KINCI
N., A
KINCI
G., T
ESCHNER
M.:Unified spray, foam and air bubbles for particle-based fluids.
The Visual Computer 28 , 6 (2012), 669–677. doi:10.1007/s00371-012-0697-9 . 1, 2[IBAT11] I
HMSEN
M., B
ADER
J., A
KINCI
G., T
ESCHNER
M.: Anima-tion of air bubbles with SPH. In
Proceedings of the International Confer-ence on Computer Graphics Theory and Applications (2011), pp. 225–234. doi:10.5220/0003322902250234 . 2[ISN ∗
20] I
SHIDA
S., S
YNAK
P., N
ARITA
F., H
ACHISUKA
T., W
OJTAN
C.: A model for soap film dynamics with evolving thickness.
ACMTransactions on Graphics 39 , 4 (2020), 31:1–31:11. doi:10.1145/3386569.3392405 . 2[KBST19] K
OSCHIER
D., B
ENDER
J., S
OLENTHALER
B., T
ESCHNER
M.: Smoothed particle hydrodynamics techniques for the physics basedsimulation of fluids and solids. In
Eurographics 2019 - Tutorials (2019),The Eurographics Association. doi:10.2312/egt.20191035 . 8[KDBB17] K
OSCHIER
D., D
EUL
C., B
RAND
M., B
ENDER
J.: Anhp-adaptive discretization algorithm for signed distance field genera-tion.
IEEE Transactions on Visualization and Computer Graphics 23 ,10 (2017), 2208–2221. doi:10.1109/TVCG.2017.2730202 . 6[KLKK12] K IM P.-R., L EE H.-Y., K IM J.-H., K IM C.-H.: Con-trolling shapes of air bubbles in a multi-phase fluid simulation.
The Visual Computer 28 , 6 (2012), 597–602. doi:10.1007/s00371-012-0696-x . 2[KTT13] K IM T., T
ESSENDORF
J., T
HUEREY
N.: Closest point turbu-lence for liquid surfaces.
ACM Transactions on Graphics 32 , 2 (2013),15:1–15:13. doi:10.1145/2451236.2451241 . 8[MBT ∗
15] M
ERCIER
O., B
EAUCHEMIN
C., T
HUEREY
N., K IM T.,N
OWROUZEZAHRAI
D.: Surface turbulence for particle-based liquidsimulations.
ACM Transactions on Graphics 34 , 6 (2015), 202:1–202:10. doi:10.1145/2816795.2818115 . 2[MLJ ∗
13] M
USETH
K., L
AIT
J., J
OHANSON
J., B
UDSBERG
J., H EN - DERSON
R., A
LDEN
M., C
UCKA
P., H
ILL
D., P
EARCE
A.: Open-VDB: An open-source data structure and toolkit for high-resolutionvolumes. In
ACM SIGGRAPH 2013 Courses (2013), pp. 19:1–19:1. doi:10.1145/2504435.2504454 . 8[MMR13] M
ACDONALD
C. B., M
ERRIMAN
B., R
UUTH
S. J.: Simplecomputation of reaction–diffusion processes on point clouds.
Proceed-ings of the National Academy of Sciences 110 , 23 (2013), 9209–9214. doi:10.1073/pnas.1221408110 . 2[MRWE20a] M
ORGENROTH
D., R
EINHARDT
S., W
EISKOPF
D.,E
BERHARDT
B.: Application and sample scenes for efficient 2D sim-ulation on moving 3D surfaces, 2020. doi:10.5281/zenodo.4009208 . 8[MRWE20b] M
ORGENROTH
D., R
EINHARDT
S., W
EISKOPF
D.,E
BERHARDT
B.: Source code for efficient 2D simulation on moving 3Dsurfaces. https://github.com/dimo3d/Cappucino , 2020. 2[Red70] R
EDLICH
O.: Intensive and extensive properties.
Journalof Chemical Education 47 , 2 (1970), 154–156. doi:10.1021/ed047p154.2 . 3[RKEW19] R
EINHARDT
S., K
RAKE
T., E
BERHARDT
B., W
EISKOPF
D.: Consistent Shepard interpolation for SPH-based fluid animation.
ACM Transaction on Graphics 38 , 6 (2019), 189:1–189:11. doi:10.1145/3355089.3356503 . 8 c (cid:13) (cid:13) orgenroth et al. / Efficient 2D Simulation on Moving 3D Surfaces [RM08] R UUTH
S. J., M
ERRIMAN
B.: A simple embedding method forsolving partial differential equations on surfaces.
Journal of Compu-tational Physics 227 , 3 (2008), 1943–1961. doi:10.1016/j.jcp.2007.10.009 . 2, 6[SABK08] S
EYCHELLES
F., A
MAROUCHENE
Y., B
ESSAFI
M., K EL - LAY
H.: Thermal convection and emergence of isolated vortices insoap bubbles.
Physical Review Letters 100 , 14 (2008), 144501. doi:10.1103/PhysRevLett.100.144501 . 8[SGGM06] S UD A., G
OVINDARAJU
N., G
AYLE
R., M
ANOCHA
D.: In-teractive 3D distance field computation using linear factorization. In
Pro-ceedings of the 2006 Symposium on Interactive 3D Graphics and Games (2006), pp. 117–124. doi:10.1145/1111411.1111432 . 6[Sta03] S
TAM
J.: Flows on surfaces of arbitrary topology.
ACM Transac-tions on Graphics 22 , 3 (2003), 724–731. doi:10.1145/1201775.882338 . 2[Sto90] S
TONE
H.: A simple derivation of the time-dependentconvective-diffusion equation for surfactant transport along a deforminginterface.
Physics of Fluids A: Fluid Dynamics 2 , 1 (1990), 111–112. doi:10.1063/1.857686 . 9[SY04] S HI L., Y U Y.: Inviscid and incompressible fluid simulation ontriangle meshes.
Computer Animation and Virtual Worlds 15 , 3-4 (2004),173–181. doi:10.1002/cav.19 . 2[TFK ∗
03] T
AKAHASHI
T., F
UJII
H., K
UNIMATSU
A., H
IWADA
K.,S
AITO
T., T
ANAKA
K., U
EKI
H.: Realistic animation of fluid withsplash and foam.
Computer Graphics Forum 22 , 3 (2003), 391–400. doi:10.1111/1467-8659.00686 . 1, 2[TSS ∗
07] T
HUEREY
N., S
ADLO
F., S
CHIRM
S., M
ÜLLER -F ISCHER
M., G
ROSS
M.: Real-time simulations of bubbles and foam withina shallow water framework. In
Proceedings of the 2007 ACMSIGGRAPH/Eurographics Symposium on Computer Animation (2007),pp. 191–198. doi:10.2312/SCA/SCA07/191-198 . 2[XB14] X U H., B
ARBI ˇC
J.: Signed distance fields for polygon soupmeshes. In
Proceedings of Graphics Interface (2014), pp. 35–41. 6[XZ03] X U J.-J., Z
HAO
H.-K.: An Eulerian formulation for solv-ing partial differential equations along a moving interface.
Journalof Scientific Computing 19 , 1 (2003), 573–594. doi:10.1023/A:1025336916176 . 2[ZB05] Z HU Y., B
RIDSON
R.: Animating sand as a fluid.
ACM Transac-tions on Graphics 24 , 3 (2005), 965–972. doi:10.1145/1073204.1073298 . 6 c (cid:13) (cid:13)(cid:13)