An explicit structure-preserving numerical scheme for EPDiff

We present a new structure-preserving numerical scheme for solving the Euler–Poincaré Differential (EPDiff) equation on arbitrary triangle meshes. Unlike existing techniques, our method solves the difﬁcult non-linear EPDiff equation by constructing energy preserving, yet fully explicit , update rules. Our approach uses standard differential operators on triangle meshes, allowing for a simple and efﬁcient implementation. Key to the structure-preserving features that our method exhibits is a novel numerical splitting scheme. Namely, we break the integration into three steps which rely on linear solves with a ﬁxed sparse matrix that is independent of the simulation and thus can be pre-factored. We test our method in the context of simulating concentrated reconnecting wavefronts on ﬂat and curved domains. In particular, EPDiff is known to generate geometrical fronts which exhibit wave-like behavior when they interact with each other. In addition, we also show that at a small additional cost, we can produce globally-supported periodic waves by using our simulated fronts with wavefronts tracking techniques. We provide quantitative graphs showing that our method exactly preserves the energy in practice. In addition, we demonstrate various interesting results including annihilation and recreation of a circular front, a wave splitting and merging when hitting an obstacle and two separate fronts propagating and bending due to the curvature of the domain.


Introduction
Structure-preserving integrators [HLW06] play an important role in many applications in pure and applied mathematics.These methods are unique in that they solve a differential equation while keeping certain geometrical invariant properties of the system.Numerical integrators that respect the underlying invariants of the dynamics are especially important in computer graphics, where gain/loss of vorticity [AWO * 14] or kinetic energy [AVW * 15, MCP * 09] may lead to inaccurate long-integration and undesirable visual artifacts.The goal of this paper is to suggest an efficient structure-preserving numerical scheme for solving the challenging EPDiff equation.
EPDiff arises in several applications and domains.For example, in computational anatomy, EPDiff are the governing equations in the diffeomorphic deformation framework [MTY02], where the goal is to compute a non-linear deformation between a given pair of source and target images.In fluid mechanics, EPDiff is used in the context of turbulence with the so-called Navier-Stokes-α model.A similar formulation was devised to better handle an effect known as energy cascading [LMH * 15].Perhaps closest to our focus is the link to shallow water wave dynamics, where EPDiff are used to model ocean wavefronts, 100-200 km in length (see the images in [HS13]).We also use EPDiff to model interacting surface waves.
Intuitively, EPDiff describes the evolution of "waves" on a fixed domain that can reconnect after collision, i.e., concentrated waves can pass through each other and maintain their shape.To capture the reconnection effect, EPDiff encodes the following key idea.Upon collision, an exchange of momentum occurs between the involved fronts.Figure 2 illustrates this behavior for 1D singular waves (peakons) where the initial peaked solitons "switch" places due to exchange of momentum, and we further show a corresponding 2D example in Figure 11.The derived model can be interpreted as the advection of the concentrated momentum over the velocity field (a non-linear term) where the velocity is a smoothened version of the momentum (a non-local term).Interestingly, the above setup is reminiscent of the vorticity equation (see e.g., [Saf92]) which governs the kinematics of ideal incompressible flows where vorticity is being transported by the velocity and these quantities are linked through the Biot-Savart law.We provide a further elaborate discussion on the similarities between the models in Section 2.
The Euler-Poincaré (EPDiff) equation has received increasing attention in the literature (see e.g., [HS13]).However, while there are some works which manage to discretize this complex PDE in various configurations, its discretization is still considered a nontrivial task since it involves several challenges.For example, any discrete method must be relatively accurate because the advected momentum is highly concentrated and thus discretization errors become visible quite quickly.Moreover, the stability and behaviour of the fronts depend heavily on the particular non-linearity of the equation; if that is not discretized correctly, the concentrated waves will be unstable regardless of how numerically accurate the discretization is.Similarly, time integration is equally important as it is expected to preserve the properties of the continuous problem as much as possible.Finally, the spatial differential operators should gracefully handle boundaries and deal with curved domains.
To better assess the proposed approach, we will try to classify it with respect to other methods for simulation of fronts and waves.From a broader perspective, our model can be considered as part of the family of shallow water equations (see e.g., [Vre13]).In these models, the assumptions of columnar motion and averaged velocity over the fluid height naturally lead to a reduction of dimensionality.Thus, although the 3D Navier-Stokes (NS) equations could capture the effects we are interested in, the involved computational cost is Figure 3: Tracking fronts shown here as red curves can be done using the eikonal eq.leading to viscosity solutions (left).Alternatively, periodic functions yield superimposed interaction between waves denoted by the blue curve (middle).Our method tracks concentrated waves which interact with other waves (right).See also Fig. 5. too prohibitive for practical uses when compared to a 2D model such as ours [Bri08].Moreover, qualitative properties are usually hard to infer from the general NS model.For instance, certain singular solutions are known to exist for our model, allowing for a better qualitative and quantitative understanding of the model.
Alternative approaches to wave simulation are commonly based on insights from linearized water wave theory (see e.g., [DD91]).At the core of the linearized model we are given descriptions of the wave function as a sum of sinusoidal functions and of the wave propagation speed as the solution of the eikonal equation.Numerical methods based on this framework can be categorized into the following two groups.The first group of techniques treat the propagated front as a geometrical wave and thus obtain viscosity solutions, i.e., concentrated fronts are possible, however, the superposition principle of waves does not hold, see Fig. 3, left.Schemes from the second group approximate the wave function, and the obtained results do exhibit superposition, but modeling concentrated waves is difficult, Fig. 3, middle.Our method enjoys the advantages of both approaches, thus achieving superimposed concentrated waves as we illustrate in Fig. 3, right.
In this work, we present a fully (time-and space-) discrete scheme for the EPDiff equation, with the explicit aim that it is structure-preserving, i.e., that it conserves a physically appropriate energy, and that it is applicable to general meshes of complex surfaces.The first goal is achieved via a novel time integrator, which despite being fully explicit, is energy preserving (and therefore stable) by construction, under some mild conditions.Furthermore, the scheme is based on standard discrete (differential and interpolation) operators that are well-behaved on unstructured triangulations of curved surfaces.Overall, our method is efficient as it involves, per step, at most three sparse linear solves with a fixed matrix, that depends only on the mesh and thus can be pre-factored for the entire simulation.Moreover, as a post-processing step, we can couple the wavefronts tracked by our scheme with periodic waves of various speeds, to achieve a combined reconnection-superpositiondispersal effect.Finally, we show several results including bending of waves due to curvature effects, plausible behavior of waves interacting with boundaries, and reconnection of concentrated fronts after collision.

Related Work
Mathematical mechanics.The EPDiff equation has received increasing attention over the last decade, and the interested reader is referred to the book [HSSE09] which presents theoretical aspects of the problem, as well as several applications.Unfortunately, there exist only a few numerical methods which discretize this equation, possibly due to its complexity and related challenges.For instance, it is known [HM05] that solutions might develop discontinuities, even though the initial conditions are smooth (we also show it in, e.g., Fig. 10).Nevertheless, various particle methods approximate the solution as a linear combination of Dirac functions and exhibit momentum preservation [CTM12] or apply different boundary conditions [CKL14].The Eulerian approach presented in [HS13] uses an explicit Runge-Kutta temporal integrator on logically rectangular flat domains, and thus, it is not structure-preserving.Probably closest to our method is the technique which was recently suggested by Larsson et al. [LMMM16].They introduce structurepreserving schemes for EPDiff on a regular structured grid.Their method preserves energy but it is not time invertible.In addition, extending their method to support curvature effects is non-trivial, as it would require generalizing it to handle triangle meshes which are frequently used for representing curved domains.
Fluid simulation.The topic of ocean and wave simulation is well researched in computer graphics and reviewing the vast associated literature is beyond the scope of this paper.In our discussion, we mainly focus on methods which specifically model waves and their interaction, and we refer to the available reviews [Bri08,DCGG11] which provide an extensive discussion.Moreover, we emphasize that we only suggest to capture wave-like phenomena using EPDiff as a testing environment for our numerical scheme.In practice, while our approach may provide certain advantages over existing techniques, it should be generally considered as a complementary approach.We base our summary of related-work on the following classification.In the linear model of waves, the function whose gradient determines the propagation speed is known as the phase function.Mathematically, describing the aforementioned reconnection effect is equivalent to requiring a multi-valued phase function, whereas single-valued phase functions discard wave interactions, i.e., wave reconnection is lost.In what follows, we divide methods depending on whether their approximated phase function is multivalued or single-valued.
A popular approach for phase function approximation is to assume it is a sum of fixed propagation velocities [MWM * 87], which can be further improved [T * 01,HNC02].A natural extension of the former approach aims for a better integration of the eikonal equation by assuming varying speeds [FR86,Pea86], and consequently, it allows for additional wave effects.Another common alternative is to employ reductions such as shallow water equations [KM90] whose extensions enable intricate effects of interaction with rigid and soft bodies [CM10] and even breaking of waves [TMFSG07].Keeler and Bridson [KB14] simulate deep ocean waves by solving for a potential flow using a boundary integral equation.Yuksel et al. [YHK07] present an interesting approach where fronts are represented using particles achieving interaction with dynamic objects at realtime rates.Other works attempt to design waves using guide shapes [NB11] or through fitting a desired look [NSB13].
Supporting multi-valued phase functions can be achieved through solving the volumetric Navier-Stokes eqs.[FF01].However, while this method captures several effects and in particular concentrated fronts are doable, the involved computational cost is too restrictive when compared to other methods.iWave [Tes04] offers an attractive framework for achieving wave behaviors as diffraction and refraction.Recently, Jeschke and Wojtan [JW15] facilitated a Lagrangian technique known as wavefront tracking to reconstruct high-frequency wave functions, and to obtain interesting wave behaviors while improving on former methods [GLS00, TB87, GM02].
Computational imaging.Finally, we briefly mention that the EPDiff equation also appears in computational anatomy in shape matching problems.Given initial and final outlines of an image, the goal is to interpolate between the input outlines using an optimization problem whose solution satisfies the EPDiff equation.An attractive application of this framework is the 3D reconstruction of the human brain from 2D PET scans [BMTY05].However, since the above setup involves boundary value problems, it requires other approaches than those used for initial value problems as ours, see e.g., [AVBC16].Therefore, we only point out this interesting link and refer to a detailed review on the subject [MTY02].

Contributions
The contributions of our method can be summarized as follows.
Discrete EPDiff on curved domains.Our method is based on an efficient discretization of the EPDiff PDE by devising a novel fully explicit splitting scheme for the temporal integration which is energy preserving by construction.Additionally, our method offers a tradeoff between exactly maintaining the velocity-momentum nonlocal relation or obtaining a time invertible scheme.Finally, our spatial discretization makes use of standard operators defined directly on triangle meshes and thus can be easily implemented.O. Azencot, O. Vantzos & M. Ben-Chen / An explicit structure-preserving numerical scheme for EPDiff Concentrated reconnecting wavefronts.We simulate geometrical fronts whose behavior closely reflects the non-linear interaction between waves.Our fronts are different from the viscosity solutions arising from the eikonal equation.In addition, we convolve the simulated data with various wave sources to support globallysupported periodic waves in a fashion similar to wavefront tracking methods.
Non-trivial effects.We show several intricate results on flat and curved domains including annihilation and recreation of a circular front (Fig. 10), a wave splitting to two as it interacts with an obstacle (Fig. 15), waves exhibiting different spatial profiles due to the underlying curvature (Fig. 12) and other interesting results.

Smooth Setting
We seek to simulate, via the EPDiff equation, the motion of wavefronts such that the resulting fronts exhibit profiles which are realistic, see e.g., Fig. 5.Although the equation itself is the product of a rather complicated chain of reductions, regularizations and approximations of the Navier-Stokes equations, we will try to motivate it here in a more direct manner.
The key intuition is that an expanding wavefront can be seen as a self-propagating concentration of linear momentum, much in the same way that a vortex is a self-propagating concentration of vorticity (i.e., angular momentum).In other words, a wavefront forces the fluid to move linearly, just as the vortex forces the fluid to rotate around it.An initial set of vortices, each localized along a curve (vortex filament) or a hypersurface (vortex sheet), interact in a nonlocal manner as they are advected by the sum of their individual induced velocity fields.In the case of an inviscid and incompressible fluid (the Euler equation regime), this motion is described by the vorticity equation (see e.g., [ZBG15]), which can be written as This PDE is used to update the vorticity, given a velocity field, via a combination of transporting and stretching.To close the system we need to specify the relation between the velocity and the vorticity, which is given by ω = ∇ × v.In integral form, this relation is the Biot-Savart law v(x) = 1 4π ω(y)×|x−y| |x−y| 3 dy, which gives the velocity as a function of the vorticity, and thus makes clear the non-local interaction between the vortices.
In a similar manner, the EPDiff equation [HSSE09, Eq. 11.20] is given by together with the non-penetration boundary condition v • n = 0 at the boundary of the domain, and it describes the motion of wavefronts, i.e., concentrations of linear momentum localized along hypersurfaces, each with its own direction and speed.Like the vorticity, the momentum m is advected by the velocity field v, while the stretching terms (with an extra (div v)m term because of compressibility) are necessary to ensure the preservation of the (total) kinetic energy 1 2 v • m dx.As before, we need to close the system with a suitable constitutive relation between the velocity v and the momentum m.In most applications, where the fluid density ρ is constant, we essentially identify the momentum with the velocity; the Navier-Stokes for instance, although describing momentum transport as EPDiff does, are written exclusively in terms of the velocity.However, given the almost singular distribution of momentum in this case, attempting to self-advect with v = m does not lead to a well-behaved flow.We take instead the Kelvin-filtered velocity (see e.g., [FHT01]): where ∆ is the vector Laplacian and α is a regularization parameter that regulates the typical thickness of the fronts.The non-local nature of the velocity field, and hence of the interaction between the wavefronts, becomes clearer if we rewrite it in the integral form v(x) = G(x, y)m(y) dy where G(x, y) is the Green function of the elliptic operator Dα.In Fig. 6, we demonstrate the non-local relationship between the velocity and the momentum, by showing how initial conditions are typically set for varying values of α.
Our goal is to design a stable and efficient numerical scheme for EPDiff, and our guiding principle is to use an as-cheap-as-possible integrator which still satisfies the properties of the smooth equations: energy preservation, time-reversibility and the constitutive relationship between the momentum and the velocity.To do that, we re-formulate EPDiff in a way that makes it easier to see how the these properties arise in the smooth case.Then, we will use these insights to guide our choice of numerics.
First we note that Equation (2) can be written as: Figure 5: Our method propagates momentum over velocity.We show the norm of the velocity v using the colorbar on the right.The obtained fronts exhibit a unique behavior.In particular, after collision, we obtain both the standard viscosity solution and additional fronts which appear in the periodic solution.See also Fig. 3.  with . This formulation allows us to focus on the two important operators R and Dα, and the properties they need to have for the smooth properties to hold.Specifically, we have that: Lemma.Let (m(t), v(t)) be a solution of Equation (4).If Dα is a self-adjoint operator with respect to the L2 inner product u, v := u • v dx, i.e., u, Dαv = Dαu, v , then d dt 1 2 v, m = 0. Proof Since Dα is time-independent and self-adjoint, we get mt = Dαvt and vt , m = v, mt where mt is the time derivative of m, and so it follows that d dt Using Eq. ( 4), we have v, mt = v, R(v, m) .Thus to complete the proof we need v, R(v, m) = 0 for any v, which can be shown using standard vector calculus identities (see Appendix A).
Thus, we have identified the two properties which are necessary for energy preservation, namely Dα should be self adjoint and v, R(v, m) = 0 for any v.We will therefore assume these properties as given, and design a time discretization accordingly.We first propose a structure-preserving integrator (SPI), solving separately for v and m, such that the resulting scheme is time-invertible at the price of allowing v and m to drift from their constitutive relation.We then show how this can be amended by giving up time-invertibility and proposing an averaged structure-preserving integrator (ASPI).We also mention how all these properties can be achieved with an implicit scheme, at the expense of a higher computational burden (see Appendix B).

Time discretization
Structure-preserving integrator (SPI).We integrate the abstract form of the EPDiff equation as given in Eq. (4) using a splitting scheme.Our integrator is made up of successive substeps that update either the momentum or the velocity alone, and it eventually outputs the next (m k+1 , v k+1 ) at time t k+1 = t k + τ.While the proposed method preserves the energy and can be inverted in time, it unfortunately breaks the relation between m and v as is defined in Eq. (3).Thus, we introduce the notation v = D −1 α m and m = Dαv to simplify notation and to better distinguish between the quantities m, v that we track in practice vs. their corresponding v, m.Our update rule is then where (m k , v k ) are the momentum and velocity at time t k .It is straightforward to check that the proposed scheme is fully explicit and second order accurate in time.The proposed form of "half step-entire step-half step" can be identified on the one hand as a Strang splitting [Str68], while on the other hand it closely mirrors the structure of the well-known second order symplectic Verlet integrator [Rut83] for separable Hamiltonian systems.
Energy preservation and time invertibility.We can show that each substep, and therefore the entire scheme, leaves the kinetic energy E = 1 2 v, m invariant, see the Appendix A for a proof.In Fig. 7 we show that indeed, energy is preserved in all our simulations up to machine precision.A second important property which Eq. (2) satisfies is time invertibility, namely that the transformation m → −m, v → −v and t → −t leaves the constitutive relation (3) and the entire PDE invariant.This property also holds in the time discrete case when integrating with the update rule in Eqs.(5) in the following sense.Given a set of solutions {(m i , v i )} computed using SPI, one can arrive at (−m k , −v k ) by starting from (−m k+1 , −v k+1 ) and applying the integrator.This can be shown for a substep and similarly for all steps by noticing that −m k+ 1 2 = −m k+1 + τ 2 R(−v k+1 , − mk+1 ) since R is bilinear and Dα is linear.See Fig. 8 for an example of this property.Second order convergence.Our splitting scheme voids the connection between velocities and momenta, i.e., v k+1 is not exactly the smoothened version of m k+1 and a certain drift occurs.Nevertheless, we confirm that our method is second order in time by measuring the drift and showing convergence in Fig. 9.We tried two different and complex scenarios on the sphere (red and blue) with exponentially decreasing time steps τ.We quantify the numerical error using the α norm u 2 α := u, Dαu : that is fully implicit.However, this scheme is less practical since it involves a linear solve per step of a matrix which depends on the current step (compare with SPI which only inverts Dα).Instead, we show next a modified integrator that is again fully explicit, and offers a practical tradeoff -the energy and the constitutive relation m = Dαv are preserved, but time invertibility is no longer maintained.
Averaged structured-preserving integrator (ASPI).Using the integrator (5) we have at the end of a full step v k+1 = D −1 α m k+1 .One could attempt to rectify this by, for instance, discarding the computed velocity and taking the associated quantity vk+1 := D −1 α m k+1 instead.The problem is that, although this restores the constitutive relation, it breaks energy preservation since Interestingly, the tuple ( vk+1 , m k+1 ) overshoots the energy, whereas the alternative tuple (v k+1 , mk+1 ) undershoots it (see the section 'Energy preservation for the ASPI' in Appendix A).Thus, since v k+1 and vk+1 are both estimates of the velocity at time t k+1 , we can look for a (step dependent) convex combination that blends these velocities and keeps the energy fixed.We therefore choose a scalar parameter λ ≡ λ( vk+1 , v k+1 ) as and compute the next velocity using v λ = (1 − λ)v k+1 + λ vk+1 , with its associated momentum m λ := Dαv λ .It turns out (see Appendix A) that when the ratio of the inner products is positive, then v λ , m λ = v k , m k and the energy is preserved exactly.This is indeed the case under normal circumstances, since v and v are both approximations of v(t k+1 ).If, on the other hand, the ratio is negative, we still have v λ , m λ = v k , m k + O(τ 3 ).

Spatial discretization
Geometry, functions, vector fields and inner products.We consider a triangle mesh M where V, E and F are its vertices, edges and faces, respectively.When required, we attach subscripts to  quantities to denote their domain, e.g., f V is a function on vertices, and we similarly use E or F. We use a typical setup of piecewiselinear functions and piecewise-constant vector fields, along with their associated inner products.Namely, we represent real-valued functions as scalars on the vertices of the mesh, i.e., f ∈ R |V| , and extend them to the whole mesh using piecewise-linear hat basis functions.Notice that in this case gradients of functions are piecewise-constant, and thus in our setup vector fields are sampled using a constant tangent vector per face, i.e., v ∈ R 3|F | .
For defining mass matrices we require vertex, edge, and face areas, denoted by a V ∈ R |V| , a E ∈ R |E| and a F ∈ R |F | , respectively.Vertex areas are 1/3 of the total area of their adjacent triangles, and similarly, we use 1/3 of the sum of area of the neighboring triangles for edge area.The resulting diagonal mass matrices are denoted by for vertices, edges and faces, respectively.Additionally, we use , to denote L 2 inner products; for example the inner product between functions on vertices is given by f For L 2 inner products between vector fields on the faces and on the vertices, we define G 3F , G 3V by repeating the corresponding mass matrices 3 times.This yields, e.g., the discrete inner product v, m = v T G 3F m.
We define a matrix I F V ∈ R |V|×|F | which interpolates quantities from faces to vertices, i.e., I F V (i, j) = (1/3)a F ( j)/a V (i), iff vertex i belongs to face j and 0 otherwise.Similarly, I V F interpolates data from vertices to faces and is defined by To reduce clutter we sometimes denote an interpolated quantity Differential Operators.To preserve energy discretely, we need the operators grad ∈ R 3|F |×|V| and div ∈ R |V|×3|F | to fulfill discrete integration by parts, thus we define div = −G −1 V grad T G 3F .This definition coincides with the standard construction of these operators (e.g., as defined in [BKP * 10, Chapter 3]).To construct Dα we further need a vectorial Laplacian.Instead of using the scalar Laplacian component-wise, we use the following Hodge Laplacian [dGDT15] suited for curved domains where J is an operator which rotates tangent vectors with respect to the face normal, and div E is the edge-based divergence operator.
It is easy to verify that ∆ is self-adjoint with respect to the inner product since the operator in the parentheses can be written as a multiplication between a matrix and its transpose.
Finally, we need an operator ∇ m which provides the Jacobian of a vertex-based vector field m, at the j-th face.Thus we define (∇ m) j = (grad mx) j (grad my) j (grad mz) j T ∈ R 3×3 .
Discrete EPDiff equation.Using the above notation and discrete operators, we can discretize R(v, m) as follows: Comparing with the EPDiff equation (2), this is a rather direct discretization, additionally using necessary interpolations between face-and vertex-based quantities (see also the section 'Energy Preservation for the SPI' in Appendix A).As we have seen in the previous sections, the kinetic energy is preserved as long as v, R(v, m) = 0 for all (discrete) velocities v. Indeed, we note that v, R(v, m) = ∑ j a F ( j)v T j R(v, m) j and the contributions of the first two terms in (7 Likewise, the contribution of the last two terms also cancel out.In the smooth setting, integration by parts leads to v, grad(m • v) = − v, (div v)m .Similarly, using our discrete operators, it is straightforward to show that v, grad( m • ṽ) = − v, I V F (div v) m .Furthermore, the resulting scheme is timeinvertible, as it is easy to check that R(−v, −m) = R(v, m).

Implementation details
We implemented our technique in MATLAB.Given initial m 0 and its corresponding v 0 , at each step we integrate in time using our SPI or ASPI methods by computing the discrete R(v, m) which appears in Eq (7), evaluating the discrete Dα and facilitating its precomputed factorization for solving the required linear systems.
The initial momentum is set by prescribing a delta function over the faces which represent the front, and smoothening it using a small parameter β.This provides the norm of m 0 and its direction was commonly chosen as the normalized gradient of a distance function.Then, the initial velocity v 0 is computed by the non-local relation, see also Fig. 6.A rule of thumb for choosing α is to use 3-4 times the average edge length of the mesh.We provide the parameters α and β for our simulations in Table 1.
Per step, the inverse of the non-local Dα matrix is needed twice in SPI and three times in ASPI.Fortunately, this matrix depends only on the mesh and is independent from the entire simulation.Thus, we pre-factor it once as a pre-processing step using the SuiteSparse library [Dav06].With the pre-factorization in place, the estimate for the computational cost per step can be roughly divided to 60% for evaluations of R and 40% for the linear solves for SPI and 50%/50% for ASPI.We show further details and timings of our simulations in Table 1.
Finally, we employ a basic Courant-Friedrichs-Lewy (CFL) rule for the dynamic time stepping.Specifically, the next time step τ k+1 is updated using the relation where γ is a fixed constant in all of our simulations taken to be 0.4 for SPI and 0.8 for ASPI, and the average edge length δx is divided by the infinity norm of the current velocity, so that propagation distance is limited to at most a single cell per step.
Limitations.Our method has a few inherent disadvantages which suggest future exploration.For instance, our model makes use of the averaged velocity per fluid column, thus effects that result from different velocities, e.g., breaking of waves, are not possible.Another shortcoming of our method is its Eulerian nature, enforcing the usage of highly resolved triangle meshes with relatively small time steps.However, since the simulated fronts are extremely concentrated, one might consider to extend our method to incorporate adaptive techniques and we leave it for future investigation.

Results
Self-collapsing circular front.The initial circular front is propagating quickly towards the center of the bowl till it collapses into a small region (Fig 10, top row).The visible amount of fluid is preserved in this experiment exactly due to energy conservation which is crucial in such stress tests.However, the velocity-momentum relation must be also maintained in such extreme cases.For comparison, the SPI method stops (τ becomes zero) at the collapsing point, whereas the ASPI method passes this point, allowing the circular front to re-emerge and flow outward (Fig 10, bottom row).
Newton's cradle on a torus.In Fig. 11, three rings of different momenta are placed on the torus where the ring with the higher momentum travels faster than the other two rings.When two rings collide, an exchange of momentum occurs, similar to the 1D case shown in Fig. 2. Thus, this experiment mimics the rigid body scenario where there is a perfect exchange of momentum between swinging spheres.Notice that while momentum is not perfectly exchanged in our simulation, the general shape and intensity of the rings are preserved.
Multiple fronts on a curved geometry.We can easily set multiple fronts in the same simulation allowing them to interact between themselves.Additionally, we show that the left front exhibits a different profile compared to the right front due to the underlying bump (Fig. 12, top row, left).Then, the fronts collide and secondary waves naturally emerge (Fig. 12, top row, middle and right).Later, the concentrated waves continue to propagate and interact with the scene in a plausible manner (12, bottom row).
Flow behind obstacles.

Wavefront-driven globally-supported periodic waves
In this section we show how our approach can be easily extended to support periodic wave functions in a post-processing step.To this end, we employ a basic version of a wavefront tracking method (see e.g., [JW15]) where our simulated concentrated wave replaces the traditional tracked front resulting from solving the eikonal eq.Specifically, we stack the interpolated velocity norms onto a matrix V ∈ R |V|×|T | , i.e., each column is given by V j = I F V |v j |, where |v j | ∈ R |F | is the point-wise norm of the velocity at time t j ∈ T , where T = {0, . . .,tmax} is the set of discrete times that we have data for.Interestingly, the non-zero entries of the matrix V exactly identify the locations of our concentrated wave over time.Thus, each non-zero V i j represents that the front visited node x i at time t j .Notice that due to the nature of the EPDiff eq., the concentrated fronts are attenuated over time, therefore, in practice, we modify the entries by applying a gamma correction, (Vγ) i j = Vmax(V i j /Vmax) γ , γ < 1.
What is left to determine is the correct amplitude per node at each time t.We assume that the source of the wavefront starts oscillating at time t = 0 with an amplitude that is a sum of N sinusoidal modes with amplitudes a k and (angular) frequencies ω k .The amplitude then at node x i sampled at a time t, can be written as a convolution of a) the velocities V i j at the node at various times t j with b) the time-delayed amplitude of the source φ(t − t j ).As mentioned, every non-zero velocity entry V i j implies that a new wavefront arrives at node x i at time t j , bringing with it a time-shifted (and potentially attenuated) copy of the source signal.Taking into account the extra effect of (shallow/deep) water dispersion, where different frequencies propagate at different speeds [JW15], we arrive at the following sum over modes φ k and times t j :   Table 1: We show various statistics for our simulations including the number of vertices |V|, the average edge length δx, the smoothening parameters α, β, the number of steps and computation times in seconds.In practical terms, the calculation of the vector η(t) of the amplitude of all the nodes x i at a given time t can be performed efficiently as the matrix-matrix-vector product: where Vγ is the gamma-corrected wavefront velocity matrix, the columns of Φ(t) ∈ R |T |×N are the time-shifted sinusoidal modes sampled at the times T of the form sin(ω k max(b k t − T , 0)), and a = (a 1 . . .a N ) T is the vector of the amplitudes of the various source modes.Note that the generated wave η(t) is a piecewiselinear function on vertices which is sufficiently sampled in time and space due to our CFL condition (see Section 4), thus linear interpolation in space produces reasonable results.In Fig. 14 we show how a single circular front (left) can be convolved with a single periodic source to produce a periodic behavior (middle).Also, coupling with two sources traveling at different speeds, allows to achieve dispersion-like effects (right).

Conclusion
We proposed a novel numerical scheme for the EPDiff equation while maintaining geometrical invariants.We tested our approach in the context of simulating concentrated and globally-supported waves on general triangle meshes.Moreover, we presented two numerical splitting schemes which are fully explicit and conserve the energy by construction, allowing to choose time invertibility over preservation of the velocity-momentum relation and vice versa.Our method is extremely efficient since it involves at most three inversions of a fixed matrix which we factorize as a pre-processing step.Additionally, we showed that our numerical data can be used in the context of wavefront tracking methods where several linear wave sources are blended to obtain various wave effects such as periodicity and dispersion.We demonstrated interesting wave-wave and wave-scene interaction effects on flat and curved domains.
We believe that many potential extensions are worth exploring.In particular, adding various forcing terms to the governing PDE will probably yield intricate new effects.Also, extending the technique to support adaptive meshes will be useful, since large parts of the domain do not contain any parts of the resulting front, and thus the resolution of these parts can be reduced.Finally, allowing layers of velocities instead of a single tangential vector might be an interesting direction as it will allow to achieve a richer span of motion effects such as breaking of waves.

Figure 2 :
Figure 2: Momentum exchange.The 1D simulation illustrates how two peakons with different momenta interact over space (horizontal axis) and time (vertical axis), left.Notice that the left peakon transfers momentum to the right peakon when they collide.Paths of the left/right soliton peaks are shown with red/blue curves on the right.
c 2018 The Author(s) Computer Graphics Forum c 2018 The Eurographics Association and John Wiley & Sons Ltd.

Figure 4 :
Figure 4: Our method is adapted to work on curved geometries, allowing to propagate general initial wave profiles.
c 2018 The Author(s) Computer Graphics Forum c 2018 The Eurographics Association and John Wiley & Sons Ltd.

c
2018 The Author(s) Computer Graphics Forum c 2018 The Eurographics Association and John Wiley & Sons Ltd.

Figure 6 :
Figure6: Typical initial conditions.(left) A concentrated unit length momentum in black is smoothened with a vector Laplacian with parameter β yielding m 0 (left, middle).Then, v 0 is achieved by smoothening m 0 with parameter α (right, middle).Taking twice as large parameter leads to a smoother result (right).Notice that the direction of m 0 and v 0 is not affected by the smoothening, and it is shown with the blue arrow (left).

Figure 7 :
Figure 7: The relative energy E k /E 1 for all the simulations we demonstrate.Notice that the error is close to machine precision since ε = 2 × 10 −13 .

Figure 8 :
Figure 8: Our SPI scheme is time invertible which allows to flow forward in time (top row and left bottom row), and then flow backwards and arrive at the initial conditions (black frame to the right).

Figure 9 :
Figure 9: Quadratic convergence in time.Error plots for two simulations (blue, red) where we measure how far v is from the smoothened m for decreasing time steps.For comparison, we show a quadratic function of τ (black).See the text for details.

Figure 10 :
Figure 10: A circular front is shrinking towards the center of the bowl till it collapses, then the front re-emerges and propagates outward.
Fig. 15, top left shows how a concentrated wave breaks when it hits the cylinder.The front then passes the obstacle and reconnects afterwards, see Fig. 15, top right.It continues to sweep through the pool in a natural way, while additional fronts are generated from behind the cylinder, see Fig. 15, bottom.
The parameters b k = b(ω k ) implement the dependency of the speed of the sinusoidal waves (as a fraction of the fixed pre-calculated speed of the wavefront) on the frequency.In Fig.13we demonstrate a graphical interpretation of the proposed method, where a 1D wave is traveling and bouncing off the boundary of the domain over time (black line).The front data is then coupled with two wave sources propagating at different speeds and different frequencies (blue and red curves shown on top).The resulting convolution (black vertical curve to the right) is the sum of the separate convolutions (red and blue vertical curves).

Figure 11 :
Figure 11: Newton's cradle on a torus with rings of different momenta.Upon collision, an exchange of momentum occurs allowing the rings to "switch" places while roughly maintaining their original shape and intensity.Fig 2 shows a similar experiment in 1D.

Figure 12 :Figure 13 :
Figure 12: Multiple fronts on a curved domain propagate with different profiles due to the underlying geometry (top row, left).When the fronts meet, additional concentrated waves appear (top row, middle and right).The general behavior is realistic showing wave-wave and wave-scene interactions (bottom row).

Figure 14 :
Figure 14: A circular front (left) is coupled with a single shifted source, producing a periodic wave (middle).Using the same technique, several sources can be blended, achieving dispersion effects.