A Linear and Angular Momentum Conserving Hybrid Particle/Grid Iteration for Volumetric Elastic Contact

We present a momentum conserving hybrid particle/grid iteration for resolving volumetric elastic collision. Our hybrid method uses implicit time stepping with a Lagrangian finite element discretization of the volumetric elastic material together with impulse-based collision-correcting momentum updates designed to exactly conserve linear and angular momentum. We use a two-step process for collisions: first we use a novel grid-based approach that leverages the favorable collision resolution properties of Particle-In-Cell (PIC) techniques, then we finalize with a classical collision impulse strategy utilizing continuous collision detection. Our PIC approach uses Affine-Particle-In-Cell momentum transfers as collision preventing impulses together with novel perfectly momentum conserving boundary resampling and downsampling operators that prevent artifacts in portions of the boundary where the grid resolution is of disparate resolution. We combine this with a momentum conserving augury iteration to remove numerical cohesion and model sliding friction. Our collision strategy has the same continuous collision detection as traditional approaches, however our hybrid particle/grid iteration drastically reduces the number of iterations required. Lastly, we develop a novel symmetric positive semi-definite Rayleigh damping model that increases the convexity of the nonlinear systems associated with implicit time stepping. We demonstrate the robustness and efficiency of our approach in a number of collision intensive examples.


Impulse Iterations Convergence
We present a momentum conserving hybrid particle/grid iteration for resolving volumetric elastic collision.Our hybrid method uses implicit time stepping with a Lagrangian finite element discretization of the volumetric elastic material together with impulse-based collision-correcting momentum updates designed to exactly conserve linear and angular momentum.We use a two-step process for collisions: first we use a novel gridbased approach that leverages the favorable collision resolution properties of Particle-In-Cell (PIC) techniques, then we finalize with a classical collision impulse strategy utilizing continuous collision detection.Our PIC approach uses Affine-Particle-In-Cell momentum transfers as collision preventing impulses together with novel perfectly momentum conserving boundary resampling and downsampling operators that prevent artifacts in portions of the boundary where the grid resolution is of disparate resolution.We combine this with a momentum conserving augury iteration to remove numerical cohesion and model sliding friction.
Our collision strategy has the same continuous collision detection as traditional approaches, however our hybrid particle/grid iteration drastically reduces the number of iterations required.Lastly, we develop a novel symmetric positive semi-definite Rayleigh damping model that increases the convexity of the nonlinear

INTRODUCTION
The deformation of volumetric elastic solids is a fundamental aspect of computer graphics and related disciplines.Inertia-driven elastic deformation (particularly in response to collision and contact) adds rich detail and realism to animations.This is particularly true for animating characters with biomechanical soft tissues like muscle, fat, skin, etc. [Brunel et al. 2021;McAdams et al. 2011;Milne et al. 2016;Smith et al. 2018].Indeed, effectively conveying the transition from squash to stretch is one of the key principles of animation [Thomas and Johnston 1995].Since the pioneering work of Terzopoulos et al. [1987], researchers have developed approaches to generate elastic deformations with numerical simulation and the governing physics of elasticity.Graphics researchers often borrow and adapt techniques from the computational mechanics literature in this regard, however, despite many decades of research by both communities, the simulation of large-strain elastic solids with contact and collision remains an active area of research.
The governing physics of large strain elastic solids are comprised primarily of two factors: stress-based internal forces arising to resist deformation (often from distortion energies) [Gonzalez and Stuart 2008] and contact forces at the interface between colliding bodies [Belytschko et al. 2013; Kikuchi and Oden 1988].In both the graphics and broader computational mechanics literature, the numerical treatment of stress-based forces is relatively settled.The finite element method (FEM) discretization of spatial stress gradients is used by the vast majority of researchers due to its support of complex geometries and its generally favorable numerical properties (e.g.symmetric discretization, spectral reproduction, etc.) [Belytschko et al. 2013;Hughes 2000].In contrast, the numerical treatment of contact and collision is far less settled.Contact forces happen at such fast time scales that they are effectively discontinuous [Larionov et al. 2021;Li et al. 2020a] and as a consequence, their numerical discretization is far more delicate.
There are many existing options for resolving these terms and each of them has its relative pros and cons.Penalties, barriers, and repulsive forces effectively regularize the collision response, often idealizing it in terms of potential energy increasing with material overlap [Baraff and Witkin 1998;Barbič and James 2007;Gast et al. 2015; Moore and Wilhelms 1988;Spillmann et al. 2007;Teng et al. 2014;Teran et al. 2005].Iterative discrete contact resolution between mesh facets like points, triangles, and edges has a long history in computer graphics and can provide strong collision-free assurances [Bridson et al. 2002;Cundall and Strack 1979;Harmon et al. 2009;Müller et al. 2007;Provot 1997;Volino et al. 1995;Wu et al. 2020].Methods that formulate the problem in terms of constrained optimization are also very powerful, often resulting in Linear Complementary Problems (LCP) [Baraff and Witkin 1992;Larionov et al. 2021;Li et al. 2021;Otaduy et al. 2009;Wriggers and Laursen 2006;Wriggers et al. 1990].Hybrid Lagrangian/Eulerian methods that make use of the natural collision resolution induced by numerical discretization in an Eulerian view are increasingly popular [Fan et al. 2013;Han et al. 2019;Jiang et al. 2017a;Levin et al. 2011;McAdams et al. 2009;Yue et al. 2018 We develop an approach that uses a hybrid Lagrangian/Eulerian method combined with an iterative discrete contact model to attain the positive benefits of both approaches while minimizing their respective drawbacks.Particle-In-Cell (PIC) is a hybrid Lagrangian/Eulerian (particle/grid) technique originally developed for compressible flow applications [Harlow 1964].Recently, generalizations of PIC techniques have been used to resolve collisions with a diverse range of materials in graphics applications.Particles in a discrete deformable object on a collision trajectory have velocities that become discontinuous at the moment of collision, however, the grid transfers in PIC techniques have a regularizing effect that prevents collision by preventing these discontinuities (see [Marquez et al. 2023]).However, as noted by Han et al. [2019] and Fei et al. [2021], these techniques suffer from numerical cohesion and friction by the same mechanism.We define a conservative augury Affine-Particle-In-Cell (APIC) [Jiang et al. 2015;Tupek et al. 2021] technique that does not suffer from numerical cohesion or friction and that unlike Han et al. [2019] is perfectly linear and angular momentum conserving.We show that a notable benefit of this conservation is increased stability with large-time steps in practical simulations.PIC techniques generally require particle sampling to be comparable to grid resolution.Disregarding this constraint degrades the quality of their ability to prevent collisions (see [Marquez et al. 2023]).When viewing these transfers as operators on the boundary of a deformable object mesh, as we do, this requirement can be unrealistic.For example, when the unstructured mesh representing the deformable object has varying particle density on its boundary, the coarsest portion will determine acceptable grid resolution.This is undesirable since an overly coarse background grid leads to less effective collision resolution in practice.We design novel conservative resampling and subsequent downsampling operators for the boundary of the mesh to prevent these issues.
While our conservative hybrid approach is free from numerical cohesion and friction, maintaining this assurance can degrade its ability to prevent all collisions in practice.To resolve this, we augment our approach with the discrete contact model of Bridson et al. [2002].We show that our combined approach drastically reduces the number of iterations that would be required from their approach alone; furthermore, we attain their strong collision-free assurances.We incorporate this hybrid collision model into a novel predictor/corrector implicit time stepping scheme.Backward Euler time integration is the primary building block of our approach and we add a novel Rayleigh damping model to regulate its inherent numerical damping.Our Rayleigh damping model has a symmetric positive semi-definite linearization which increases the convexity of the nonlinear backward Euler systems and prevents the need for definiteness fixes [Teran et al. 2005] of finite element Hessians in practice.We summarize our technical contributions below: • A resampling technique that conserves the mass, the center of mass, the linear momentum, and the angular momentum of a collection of particles with APIC mass and momentum state.• A downsampling technique with analogous conservation properties are used to gather grid momenta associated with a resampled collection of particles into mass and momentum APIC state for the unresampled, original particles.• A momentum conserving augury iteration for removing numerical friction and cohesion artifacts in APIC transfers.• An implicit predictor/corrector time stepping scheme tailored to our collision response.
• A symmetric positive semi-definite Rayleigh damping model that increases the convexity of the hyperelastic backward Euler minimization problem.

RELATED WORK
Simulating contact between solids remains one of the more challenging aspects of finite deformation continuum mechanics simulations.We briefly discuss the most relevant works from the computer graphics and computational mechanics literature.There are many decades of research on the subject of contact and collision with deforming elastic objects.We refer the reader to reviews of the state-of-the-art-in computer graphics [Ascher et al. 2021;Nealen et al. 2006] as well as the course notes in [Andrews et al. 2022;Kim and Elberle 2022].In the engineering literature, early important contributions to Lagrangian finite element contact algorithms include [Belytschko and Neal 1991;Campos et al. 1982;Martins and Oden 1983;Simo and Laursen 1992;Wriggers et al. 1990].Mortar contact methods are now often preferred for their accuracy and stability [Popp et al. 2010;Puso and Laursen 2004], but can have high computational expense due to additional Lagrange multiplier unknowns, typically do not exactly conserve angular momentum, and not all contact collisions are guaranteed to be detected.Generally, techniques for resolving collisions with deformable bodies fall into three categories: penalties/repulsive forces [Baraff and Witkin 1998;Barbič and James 2007;Gast et al. 2015;Moore and Wilhelms 1988;Spillmann et al. 2007;Teng et al. 2014;Teran et al. 2005], iterative discrete contact [Bridson et al. 2002;Cundall and Strack 1979;Harmon et al. 2009;Müller et al. 2007;Provot 1997;Volino et al. 1995;Wu et al. 2020] and constrained optimization [Baraff and Witkin 1992;Larionov et al. 2021;Li et al. 2021;Otaduy et al. 2009;Wriggers and Laursen 2006;Wriggers et al. 1990].However, hybrid approaches are also possible [Daviet 2020;Fan et al. 2013;Han et al. 2019;Jiang et al. 2017a;Levin et al. 2011;McAdams et al. 2009;Yue et al. 2018].Our approach utilizes a combination of hybrid PIC [Harlow 1964;Jiang et al. 2015] and discrete contact [Bridson et al. 2002] models.We discuss techniques most relevant to these.Our approach is inspired by the Material Point Method (MPM), which is an extension of PIC.MPM was first proposed in [Sulsky et al. 1994] and [Sulsky et al. 1995], with the first strategy for contact introduced shortly thereafter in [Bardenhagen et al. 2000].Various improvements to MPM contact can also be found in [Homel and Herbold 2017;Huang et al. 2011;Nairn and Guo 2005;Xiao et al. 2021].A strategy for introducing discontinuities into the MPM was recently suggested in [Moutsanidis et al. 2019] and shows advantages for sliding and separating contact.Hegemann et al. [2013] make use of PIC techniques in graphics applications of ductile fracture.Here we utilize the affine particle-in-cell (APIC) method from [Jiang et al. 2015[Jiang et al. , 2017b]].The APIC augury iterations utilized in this work were first suggested in [Tupek et al. 2021], with limited examples.The initial version proved to be unstable in some situations and was unable to accurately account for large size discrepancies between foreground and background meshes.Our method successfully addresses these issues.Our approach is directly related to a few recently proposed hybrid PIC/discrete geometric collision  [Cundall and Strack 1979] with a continuum model [Klar et al. 2016;Sulsky et al. 1994] for granular materials.

BACKGROUND AND NOTATION
We represent volumetric deformable objects as a collection of particles x ∈ R  and a simplex mesh M ∈ N   × (+1) connecting them together.Our approach is designed for volumetric simulation in 3D, however, we use 2D examples for illustration.We use  = 2, 3 in our exposition to represent the case of 2D or 3D simulations respectively.We denote the boundary mesh of the volumetric object as M  ∈ N    × .Here  refers to the number of particles,   refers to the number of elements in the volumetric mesh (tetrahedra for  = 3, triangles for  = 2) and    refers to the number of elements in the boundary mesh (triangles for  = 3, segments for  = 2).We also store a velocity v ∈ R  for each particle in the volumetric mesh.Our approach requires an APIC state over the particles in the boundary of the volumetric mesh.We use A ∈ R  2   to denote the vector of affine velocities associated with particles on the boundary of the volumetric mesh.We use   to denote the number of particles on the boundary.Next, we describe details of our method related to the spatially discretized governing physics of Lagrangian hyperelasticity and backward Euler integration of their temporal dynamics.These aspects of our approach are standard in the graphics and mechanics literature and we refer the reader to Bonet and Wood [2008] and Sifakis and Barbič [2012] for more detail.The contents of the following subsections are primarily intended for a brief review and establishing notation.We use Table 3 for a quick reference to the various symbols used in the exposition.

Elasticity
We discretize and solve the partial differential equations (PDEs) of motion for volumetric hyperelastic materials [Gonzalez and Stuart 2008] with Rayleigh damping [Belytschko et al. 2013] and frictional collision constraints [Bridson et al. 2002].We use the finite element method (FEM) with linear interpolation functions over simplex meshes to integrate the PDEs [Hughes 2000].This converts spatial terms in the PDE to discrete elastic forces f  : R  → R  .These forces are the gradient of the discrete potential energy of the system PE : x ∈ R  → R. The potential energy increases as the mesh is deformed from a reference configuration (defined in terms of reference positions X ∈ R  ) to any other configuration x ∈ R  .It is defined in terms of a hyperelastic strain energy density function Ψ : R  × → R which increases with non-rigidity of the unique affine transformation that relates the reference configuration of a mesh element to that defined by x [Teran et al. 2003].This is typically referred to as the deformation gradient in the element F  : Here  0  is the measure of the element in the reference configuration, P(F) = Ψ F (F) is the first Piola-Kirchhoff stress [Bonet and Wood 2008], X  is the centroid of mesh element  in the reference configuration, x  is the position of particle  (which make up x) and    are the piecewise linear interpolation functions associated with the FEM space.Furthermore, f   (x) ∈ R  are the elastic forces on volumetric particles  that comprise the vector f  (x).Lastly, we note that we use the fixed-Corotated model of Stomakhin et al. [2012] for the energy density Ψ.

Mass
We define mass for each volumetric element  ∈ M from mass density  as   =  0  .We define the mass for each particle in the volumetric mesh by taking a portion of the mass from each element that it belongs to   =  ∈ I    +1 .Here I  is the one ring of mesh elements  ∈ M that contains the particle .We note that in a Lagrangian FEM discretization of the governing physics, conservation of mass implies that element and particle masses   do not change with the configuration of the mesh x.We use m ∈ R  to denote the vector of particle masses in the volumetric mesh.To facilitate our conservative boundary element resampling, we also need a notion of boundary element mass     for each particle 0 ≤   <  in boundary element  ∈ M  .We define the     to partition the   such that (where is the one ring of elements in the boundary mesh that contain the particle ) and |I  | is the number of elements in the set.

Backward Euler (𝐵𝐸)
Given volumetric mesh positions x  and velocities v  at discrete time   , we approximate the trajectories of the mesh vertices under the elastodynamics at the next time  +1 =   + Δ with implicit time integration.Here, Δ is the time step.We adopt the approach of Gast et al. [2015] and characterize this via the minimization Here x +1 is the backward Euler solution at time  +1 .M ∈ R  × is a lumped-mass diagonal matrix with entries equal to particle masses   .K  ∈ R  × is the Rayleigh damping matrix which we discuss in more detail in Section 8. B ∈ R   × and c ∈ R   express linear constraints over vertices in the volumetric mesh and   is the number of constrained vertices.g ∈ R  is gravity.The backward Euler velocity can be obtained after the minimization in Equation ( 1) from v +1 = x +1 −x  Δ .We will henceforth use the notation to denote the function that returns the backward Euler positions and velocities associated with the solution of the minimization in Equation ( 1).Lastly, we note that we solve the minimization problem in Equation ( 1) using Newton's method.We refer the reader to [Sifakis and Barbic 2012;Stomakhin et al. 2012;Teran et al. 2005] for more details.However, we use standard conjugate gradients to solve the linear systems.As a consequence of the increased convexity from our Rayleigh damping model, in practice we do not need to perform the definiteness fixes of the stress derivatives used Teran et al. [Teran et al. 2005].We discuss this in more detail in Section 8. Our method is comprised of three main components: backward Euler steps (Section 3.3), augury APIC hybrid particle/grid transfers (Section 5), and discrete geometric impulses (Section 7).While our discrete geometric impulse operation is the standard approach of Bridson et al. [2002], our augury APIC technique has many novel aspects.Of particular importance are our resampling/downsampling operators (Section 6) and our iterative resolution of numerical cohesion and friction (Section 5).

METHOD OVERVIEW
We first take a backward Euler step to compute tentative time  +1 positions and velocities ( x+1 , ṽ+1 ).This update is collision unaware so we rewind the positions to time   and process the tentative velocities ṽ+1 for collision.We do this in a two-step process.In the first step, we use our resampled augury APIC impulse operator to efficiently apply conservative and grid-resolution-independent impulses that prevent collision without numerical cohesion or friction (see Section 5).This defines a new boundary APIC state (v   , A +1 ).Note that after this process, the linear part of the APIC state (A +1 ) is finalized for the time step, while the constant part (v   ) will be further processed in the second step.In the second step, we apply the collision impulse model of Bridson et al. [2002] (see Section 7).We denote this process with the operator L  : R   × R   × N    × → R   .This defines the final collision processed boundary velocities v +1  .Lastly, once we have finalized the collision response on the boundary, we redo the backward Euler time step from time   to  +1 where the boundary is constrained to follow the linear trajectories determined by the finalized velocities in the previous two steps.This propagates the collision response to the interior, which allows for increased stability and larger time steps (see [Marquez et al. 2023]).We note that our collision processing has the continuous collision detection (CCD) assurances of Bridson et al. [2002] that linear trajectories of boundary points (based on their collision-processed velocities) will guarantee a collision-free state on the boundary over the time step.
We summarize this process in Fig 4 and Algorithm 1.Note that we make use of the matrix B  ∈ R   × whose rows are associated with boundary particles and whose columns are associated with volumetric particles.There is only one non-zero entry per row in the column associated with the boundary particle's location in the volumetric particle array.For convenience, we use the notation B  m to denote the masses associated with boundary particles, although this is a slight abuse of notation.

AUGURY ITERATION: L A
We leverage the natural collision prevention tendencies of PIC techniques to define our collision operator This operator maps the elastic object boundary affine momentum state to itself to resolve portions of the mesh on collision trajectories.It is composed of three components: APIC [Jiang et al. 2015] particle/grid mass/momentum transfers, novel resampling/downsampling strategies, and an augury iteration designed to prevent numerical cohesion and friction.We adopt the standard APIC transfers and refer the reader to Jiang et al. [2015] for details.Our resampling and downsampling strategies are designed to prevent the degradation of collision prevention abilities that arise with disparate boundary mesh and background grid resolutions (see Section 6).Here we discuss our conservative approach to removing numerical cohesion and friction inherent in the naive APIC transfers.
We first define the conservative operator • L RS which maps the affine velocity state of the boundary particles to itself via composition of the APIC particle-to-grid (P2G) and (G2P) operators (L P2G and L P2G respectively) together with the resampling and downsampling operations (L RS and L DS ) outlined in Section 6.We summarize this process as (2) Here we use (v  , A) to denote the APIC state of boundary (where v  ∈ R   and A ∈ R  2   are vectors containing the constant and linear parts of the APIC velocity for particles in the boundary) and (m  , p  ) to denote the grid mass and momentum state after an APIC particle-to-grid (P2G) transfer (where m  ∈ R  , and p  ∈ R  , are vectors containing the mass and linear momentum on grid nodes and  , is the number of grid node x i whose interpolation functions  i are non-zero on some boundary particle).Here,  i (x) are quadratic B-spline interpolating functions associated with grid node x i used to transfer particle quantities to grid quantities [Jiang et al. 2015;Stomakhin et al. 2013].Similarly, we use ( v , Â) and ( m , p ) to denote their resampled counterparts (where v ∈ R  N , Â ∈ R  2 N , m ∈ R N, , p ∈ R  N, , N and N, have meanings analogous to their down downsampled counterparts).
As shown in Figure 2, the APIC transfers naturally regularize particles on collision trajectories to prevent impact.However, they also prevent separation and relative sliding by the same mechanism.This was noted in Han et al. [2019] where cohesive and frictional responses were discarded in a simple and non-conservative per-boundary-particle manner.We adopt the approach of Tupek et al. [2021] to leverage the conservative-smoothing/low-pass-filter nature of the APIC operator to conservatively damp out regions where it is cohesive and/or overly frictional.We do this by iteratively applying the operator to the change it induced in the previous iteration, but only where it is deemed overly frictional or cohesive with the aim of damping the unwanted cohesion and friction to zero.By design, this process retains the momentum conservation properties of APIC techniques.We refer the reader to Jiang et al. [2015] for discussion of the conservative-smoothing/low-pass-filter properties of PIC.
We define the conservatively resampled augury APIC mapping L A by first initializing its output to be equal to the resampled APIC operator L RSA .We then interpret its output as defining impulsive changes in the boundary grid node momenta and check if these impulses are cohesive.If they are not deemed cohesive, we accept the normal component of the impulse.Next, we increment the APIC output by the difference between the conservatively resampled APIC operator L RSA applied to the cohesion processed impulse and itself.This increment does not change the momentum state since both L RSA and the identity operator are linear and angular momentum conserving.We repeat this process until we hit a maximum number of iterations or the impulses converge to a tolerance.
We reiterate that the low pass filter nature of L RSA and its repeated application to the normal component of the impulse is designed to dampen any cohesion or friction to zero.Technically, for a frictionless response, the tangential component of the impulse must be filtered along with the cohesion.However, we found that while this effectively limited these effects it degraded collision prevention.If we instead keep the tangential component whenever a change is deemed non-cohesive, we observed stronger collision prevention abilities, albeit at the expense of some numerical friction.
We outline this process in Algorithm 2. Note that the normal n  is used for determining numerical cohesion for each particle on the boundary via area-weighted normalization.Further, note that we denote the option of accepting some numerical friction for the benefit of improved collision prevention as   .

BOUNDARY ELEMENT RESAMPLING: L RS , L DS
As shown in [Marquez et al. 2023], when the resolution of the boundary mesh is too coarse relative to the background grid, the collision resolution abilities of PIC techniques degrade.In extreme cases, collisions can be completely missed altogether.We design a novel resampling strategy to remove this limitation as well as a novel downsampling strategy to return the boundary to its original resolution.
Our resampling strategy is designed on a per-boundary element manner.5).Care must be taken to conserve the total mass (   ),the center of mass (x  com ), linear (p   ), and angular momenta (l   ) of the boundary element in the resampling process.While the element mass and center of mass are defined in the standard way, i.e.
element momenta are defined with the APIC [Jiang et al. 2015] convention where boundary element particles x    contribute based on their associated grid momenta at grid node x i defined from an APIC transfer to the grid.Summing over the element points   defines the boundary element grid momentum distribution p  i =  −1   =0 p ,  i .The total linear p   = i p  i and angular l   = i x i − x  com × p  i momenta of the element (computed about the element center of mass) are defined from the element grid momentum distribution.
The total momentum distribution of the boundary is defined from the sum of the element-wise contributions p i =  p  i .The total mass and center of mass are defined similarly: The total grid linear p  and angular l  momenta (computed about the boundary center of mass) of the boundary are related to the element-wise counterparts by the relations Equation ( 3) shows that the total linear momentum of the boundary is equal to the sum of the element-wise linear momenta.As illustrated in Figure 6, the total angular momentum of the boundary is equal to the sum of the angular momenta of each element (computed about its center of mass) plus the sum of the angular momenta induced by the total linear momentum of the element and the relation of its center of mass to the global center of mass.Thus, choosing a resampling strategy that preserves   , x  com , p   , and l   gives the appropriate notion of local and global conservation.

Mass and Position Resampling: Partition of Unity and Linear Reproduction
We choose the resampled positions x    using random Poisson disc sampling within the element together with the original x    .We use the convention that the first  resampled points x    coincide with the original x    .Poisson disc sampling for the remaining points is chosen to prevent particle clumping while maintaining separation by no more than Δ 4 on average (Δ is the grid spacing).We require the resampled positions x    to be inside the element and we require that the resampled masses m    be positive and conserve the element mass m  = We also require conservation of the element center of mass x  com .Furthermore, in order to allow for momentum conservation in our downsampling strategy (see Section 6.2), we must ensure that both the center of mass of the first  resample positions and the center of mass of the remaining resample positions ( ≤   <    ) are equal to the element center of mass.We preserve the center of mass of the first  resampling points by simply scaling their original masses m    =      ,  ∈ (0, 1], 0 ≤   <  This associates a scaling of the total element mass    to these particles and the remaining resampled masses m    ,  ≤   <    must be chosen to partition the remaining element mass (1 − )   in such a way that their center of mass is equal to that of the element.This can be achieved by defining these masses in terms of interpolation functions N    ,  ≤   <    , defined over the remaining resample points as where the functions N    satisfy partition of unity (

and linear reproduction (
properties.We adopt the approach of Arroyo and Ortiz [2006] to create interpolation functions that satisfy these requirements for unstructured particles.This choice allows for total element mass conservation and the conservation of the center of mass from the partition of unity and linear reproduction properties (respectively) of the interpolating functions N    [Marquez et al. 2023].6.1.1Affine Velocity Resampling.We choose the resampled APIC velocity state v    and Â    in a manner that preserves the total linear (p   ) and angular (l   ) momenta of the boundary element and that preserves the state of the first  resample points x   = x    , 0 ≤   <  (those that coincide with the original boundary positions).Recall from Section 6 that the linear and angular momenta of the original points in the element are defined from the grid momentum distributions As outlined in the supplementary material [Marquez et al. 2023], we can relate the total linear and angular momenta of this distribution to the APIC velocity state as p   =  −1 , where  is the permutation tensor and the notation  : A denotes the vector b with indices   =      [Jiang et al. 2015 Also, the center of mass velocity is defined from the total linear momentum as v  com = 1 v    since ratios of the element vertex masses     and the total element mass   are equal to the barycentric weights of the center of mass angular momentum Here the notation  : Δ  refers to the second order tensor B with indices    =      .The corrective angular velocity Δ  affects the total angular momentum of the resampled element as Note that this shows that linearly interpolating the affine velocities A    perfectly conserves the angular momentum in the element resulting from those terms, but that linearly interpolating the linear velocities v    does generally change the angular momentum.We choose the corrective angular velocity Δ  to account for this difference: Δl   = (1 − ) The resampled velocities and affine velocities across all boundary elements together form the output state ( v , Â) in an analogous manner to the input state (v  , A).

Downsampling Affine Velocity: L 𝐷𝑆
The output m and p of the P2G operator L P2G are comprised of masses mi and momenta pi , respectively, over the grid nodes i.The grid node velocities are updated after the transfers as pi divided by mi .Resampled particles of the element can be seen as having a grid mass distribution Multiplying the grid masses by their corresponding updated grid velocities will provide the resampled element linear momentum distribution p i = m i ( pi / mi ).This has the effect of conservatively partitioning the grid momentum distribution into element-wise counterparts.
We design a conservative procedure for condensing p i back to affine velocity state v    , A    , 0 ≤   <  defined over the original element.Note that the resampled grid momentum distribution is defined over many grid nodes that do not (in general) affect the original unresampled element (see Figure 7).That is, p i is defined over grid nodes xi with  i ( x    ) ≠ 0 for some   with 0 ≤   <    .To conservatively define the APIC state over the unresampled original element, we create a grid momentum distribution p  i defined over grid nodes x i with  i (x    ) ≠ 0 for some   with 0 ≤   <  and then leverage the conservative nature of the APIC grid-to-particle (G2P) transfers from this grid momentum distribution to the affine state v    , A    , 0 ≤   < .We first compute the total  element linear and angular momenta from the grid distribution We add a portion of linear momentum p  to each grid node x i based on their associated interpolation function   and the mass ratio . We similarly distribute the angular momentum l  by computing the angular velocity    associated with grid nodes x  i assuming they have mass distribution ) and associated inertia tensor I  .We then add momenta to the grid nodes x  i associated with an angular momentum state    about the element center of mass The total mass of each distribution i ) is equal to the element mass    and their respective centers of mass are both equal to the element center of mass x  com .Furthermore, the total linear momentum of the downsampled grid distribution p  i is equal to p  and its total linear momentum is equal to l .
The details of these calculations may be found in [Marquez et al. 2023].The final output state (m  , p  ) thus consists of the grid masses  i =    i and momenta p i =  p  i .

DISCRETE GEOMETRIC IMPULSES: L I
Although our augury operator L A prevents most collisions, in practice the iterative removal of cohesive/overly-frictional terms can degrade some of its ability to prevent the collision.In Figure 8, segment penetration can occur on forced collisions with large deformations that cause penetration through the tangential component if the timesteps are large.We, therefore, augment it with an iterative impulse-based post-process using conventional point/triangle and edge/edge pairs in the boundary mesh.We adopt the approach of Bridson et al. [2002] and iteratively apply these impulses based on proximity at the beginning of the time step and based on pairs determined to be colliding during the time step using continuous collision detection (CCD).Furthermore, we adopt their friction model since by design, our augury APIC operator L A is frictionless.Particles in the point/triangle and edge/edge pairs are assumed to take linear trajectories determined from locations at the beginning of the time step and the linear part of their APIC velocity state.We use the approach of Brochu et al. [2012] for the CCD determinations.We consider this to be a process that operates on the linear part v ∈ R   of the APIC velocity state.We use L  : R   × R   × N    × → R   to denote this process.For a detailed explanation of how impulses are calculated, as well as proof that they are linear and angular momentum conserving, we refer the readers to [Marquez et al. 2023].

RAYLEIGH DAMPING
Our Rayleigh damping model serves two purposes.First, it damps non-rigid modes so that backward Euler does not suffer from purely numerical damping that reduces with smaller time steps.This adds predictive control of the amount of damping in simulations, independent of time step size.Second, it increases the convexity of the backward Euler system allowing for more rapid convergence.Rayleigh damping for FEM discretization of hyperelasticity is analogous to adding damping to a one-dimensional spring where the damping coefficient is proportional to the stiffness of the spring [Belytschko et al. 2013].However, in this case, the stiffness is given by the force Jacobian matrix f  x evaluated at the previous time step.As such our Rayleigh damping model is of the form where K  (x  ) is proportional to a modified stiffness matrix.The stiffness matrix, although symmetric, will generally have negative, zero, and positive eigenvalues.In Teran et al. [2005], an element-wise definiteness fix was used to compute an always symmetric semi-definite counterpart to admit the use of conjugate gradients for the solution of the linear systems in quasistatic problems.However, this requires the solution of a 3X3 eigensystem and three 2X2 eigensystems.Although this cost is negligible since these can be computed rapidly and in parallel, a simpler and satisfactory strategy for the Rayleigh damping model is to compute element strain energy Hessians at the polar decomposition of the deformation gradient Here   ≥ 0 is the damping coefficient where   < 1,   = 1 and   > 1 give behaviors analogous to underdamped, critically damped, and overdamped springs respectively.Note that in practice we found   = .1 to give a qualitatively appropriate amount of damping in our examples.F  = R  S  is the polar decomposition of the deformation gradient.We use the approximation P F (F  ) ≈ P F (R  ) because the matrix P F (R  ) is positive semi-definite yet still rotationally consistent.The semi-definiteness of P F (R  ) can be seen readily from the derivations in [Teran et al. 2005].From Teran et al. [2005], we know that F : P F (R  ) : F = U  FV  : P F (I  ) : U  FV  since R  = U  V  , so it suffices to show the positive semi-definiteness of P F (I  ) where I  ∈ R  × is the identity matrix.Also from Teran et al. [2005], we know that for isotropic materials, the matrix P F (I  ) is a block diagonal matrix with diagonal components A, B 12 , B 13 , B 23 such that where  and  are Lame parameters.Since each matrix on the diagonal is positive semi-definite, the matrix P F (R  ) is positive semi-definite.As shown in Equation ( 1), the Δ scaling before K  is one order lower than that of the potential energy .So in practice the mass matrix M and Rayleigh matrix K  dominate the backward Euler Hessian H  = M + ΔK  + Δ 2 f  x relative to the stiffness matrix f  x .We found that this made the H  symmetric positive definite in practice, without performing the definiteness fix on f  x .Note that without our symmetric semi-definite Rayleigh damping model, this is not the case, e.g. as was shown in McAdams et al. [2011].The collision of two circular meshes is simulated with only APIC collision responses (no L  impulses) to illustrate the grid-dependent collision resolution degradation and the numerical cohesion that arises without the augury iteration method from Tupek et al. [2021].A coarse background grid will experience numerical cohesion from grid transfers as shown in the left simulation in Figure 9 where the circles stick together.The cohesion is removed from the coarse grid with only five augury iterations.This is shown in the right simulation in Figure 9 where the circles separate after colliding.In Figure 10, we refine the background grid to highlight the importance of the resampling/downsampling strategy.Note that in the left simulation, we do not use the resampling/downsampling operators only the approach of [Tupek et al. 2021] to highlight the degradation in the collision response that occurs.The degradation arises from the coarsely spaced surface nodes relative to the background grid.In the bottom simulation, our resampling/downsampling strategy increases the foreground particle density to be high-enough relative to the background grid to capture the collision and prevent the spheres from overlapping.Thus, the collision is resolved on finer grids than would otherwise be possible.In Figure 12 two rotating elastic blocks collide to demonstrate the linear and angular momentum conservation from the resample/downsampling operators.Linear momentum is conserved through the collision for both timestep methods, however, angular momentum changes only when using backward Euler due its lack of an angular momentum conservation property [Jiang et al. 2017b].In Figure 8, two elastic blocks are compressed by a moving boundary condition.The limitations of the augury APIC method L  alone, without the geometric impulses L  are explored.The three simulations use: (left) L  alone with   = , (middle) L  alone with   =  , and (right) L  with   =  and L  .Without L  (left and middle), as compression increases, the corners and edges of the block penetrate the other block from the tangential direction.This is because the augury estimation of friction is less separated from the estimation of cohesion at regions of high curvature.The right simulation resolves the issue with the addition of L  .

Two Blocks Collision
As shown in Figure 13, our approach has the additional benefit of improving the convergence and overall run time of the impulse iterations L  .The initial frames of the collision only require a few iterations to converge as the collision occurs primarily in the normal direction, thus L  captured most of the collision.More L  impulse iterations are required when penetration in the tangential direction is allowed by L  , but significantly less than a method based purely on L  .The overall run time, including the overhead of L  is also reduced compared to L  alone.[2002] impulse iterations.This decreases the run time from 5319s to 4680s.We compare our approach with similar methods used in graphics applications.The incremental potential contact (IPC) approach of Li et al. [2020b] uses a barrier potential to perfectly resolve collisions and has been shown to handle many challenging contact problems with volumetric solids.Hegemann et al. [2013] use a similar PIC-style approach to resolving contact.We compare the performance of our approach with both approaches using a colliding spheres example (see Figure 14) over multiple spatial resolutions and with the same time step sizes (see Table 2).The PIC-style approach of Hegemann et al. [2013] performs faster than our method as the PIC/FLIP transfers are only done once as opposed to multiple times from our approach, however, this method does not conserve linear or angular momentum once collisions occur (see Figure 15).IPC performs similarly to ours, however, it is more costly at higher spatial resolution than our approach.Unlike Hegemann et al. [2013], IPC and our approach both conserve linear momentum in all cases.However, backward Euler time stepping does not conserve angular momentum (see e.g.[Jiang et al. 2017b]).When explicit (symplectic) Euler is used with our approach, conservation of the angular momentum is achieved.

Twisting Legs
In Figure 1, a body mesh is put through large deformation by imposing rotational motion constraints on the leg and arms.The collision that the twisting motion generates has a significant tangential sliding component that needs to be captured without penetration.The convergence when using only L  impulse iterations compared to using our combined approach is shown in Figure 1.The iterations required to converge are significantly reduced with our combined approach, as the overall runtime is reduced from 21082s to 15510s.

Moving Body
Our algorithm can be used to simulate collisions on animated mesh models.The motion of the interior points in a body tetrahedral mesh is fully prescribed, forcing the surface mesh to collide and respond.In Figure 11, the collisions on the surface that occur from the motion of the interior points are captured accurately with our method.

Roller Examples
In Figure 3, an armadillo mesh experiences large deformations under compression by having gravity and friction drive it through rotating cylinders.The rollers have constrained velocities and positions for each time step, the roller masses are larger than the armadillo to drive the motion.Self-collision in the armadillo and collisions between the roller and the armadillo are robustly captured.[2020b] and Hegemann et al. [2013].Furthermore, our approach drastically decreases the impulse iteration count required with the technique in Bridson et al. [2002] (see Figure 1).Note that our time step limitations are similar to Bridson et al. [2002], however, our APIC grid transfers require that the magnitude of position changes does not exceed one cell length of the background grid per time step.The APIC grid transfers were leveraged for their conservation of angular momentum, but other transfers with similar properties can be used with augury iterations.Note also that our conservative resampling strategy addresses grid-resolution restrictions on the FEM that would be present in Tupek et al. [2021] (see [Marquez et al. 2023]).In future works, the resampling and downsampling strategies similar to the one presented can be expanded to volumetric elements instead of contact surfaces.These techniques could be used in more general MPM and PIC techniques.Lastly, we mention that our use of MPM/APIC transfers can allow for the natural coupling of unstructured FEM techniques with general MPM simulations of various materials.

Fig. 1 .
Fig. 1. (Left) A body mesh is twisted by constraints at the arms and legs.Our method robustly and efficiently resolves the collision driven deformation.(Right) Impulse iteration requirements with our hybrid approach (red) are significantly reduced compared to geometric impulses alone (blue), decreasing the total simulation run time by 25 percent.

Fig. 3 .
Fig. 3.As shown here, our method is robust to extreme deformation.

Fig. 4 .
Fig. 4. Method overview: (1) collision-unaware backward Euler positions rewound to time   (2) collision processing for the boundary velocities (3) final backward Euler step with interior state reset to   and the boundary state at  +1 used as Dirichlet boundary conditions.

Fig. 5 .
Fig. 5. (Left) Boundary element  has particle quantities (blue) with implied background grid state (green).(Middle) Resampling of the particle quantities into the resample points (red) conserves various element properties.(Right) Resample points have an analogous grid representation that better represents the element in finer grids.

Fig. 6 .
Fig.6.The global angular momenta between the mesh and grid representations are equivalent.The details of the proof are in the supplemental technical document[Marquez et al. 2023], we illustrate the basis for the proof in this diagram.

Fig. 7 .
Fig. 7. (Left) The downsampling operator starts with the updated grid momenta p i .(Middle) Linear and angular momentum of the resample points' grid distribution is merged to the element center of mass: p  , l  .(Right) Momenta is distributed to the original boundary element grid in a conservative manner: p  i .

Fig. 8 .
Fig. 8. (Left) Resample augury iteration that constrains tangential sliding has segment-segment penetration.(Middle) Resample augury that allows tangential sliding has penetration as the material slides tangentially around the corners and edges.(Right) L RSA resolves the frictionless collision with no penetration.

Fig. 9 .
Fig. 9. (Left) Two circles collide using APIC transfers on a coarse grid where cohesion prevents separation.(Right) Using five augury iterations [Tupek et al. 2021] the cohesion is removed.

Fig. 10 .
Fig. 10.(Left) On a finer background grid the APIC transfers with augury iterations [Tupek et al. 2021] miss the collision.(Right) The resampling method manages to capture the collision.
Fig. 12. (Left) Collision of two rotating blocks.(Right) Linear and angular momentum plots for symplectic and backward Euler timesteps.
Fig. 13.(Left) Side view of two blocks being pushed into frictionless contact by elastic supports.(Right) Augury iteration reduces Bridson et al. [2002] impulse iterations.This decreases the run time from 5319s to 4680s.
Linear and Angular Momentum Conserving Hybrid Particle/Grid Iteration for Volumetric Elastic Contact ]. Proc.ACM Comput.Graph.Interact.Tech., Vol. 6, No. 3, Article 44.Publication date: August 2023.A . Han et al. [2019] develop a PIC-based contact algorithm for FEM meshes, including augmentation with Bridson et al. [2002] for thin strands.McAdams et al. [2009] and Yue et al. [2018] similarly use a hybrid PIC/MPM technique for hair and granular material simulations respectively.McAdams et al. [2009] augment their approach with Bridson et al. [2002], and Yue et al. [2018] combine a discrete element method (DEM) Proc.ACM Comput.Graph.Interact.Tech., Vol. 6, No. 3, Article 44.Publication date: August 2023.techniques For boundary element  ∈ M  with boundary positions x    , masses     , velocities v    and affine velocities A    where 0 ≤   < , we resample to create positions x    , masses m    , velocities v    and affine velocities Â    for 0 ≤   <    .For each element local particle index   , there is a corresponding global boundary index ; the velocity v    and affine velocity A    are then equal to v  and A  respectively, where these are the vectors which constitute the input state (v  , ) to L RS .Here we use    to denote the number of resample points in boundary element  (see Figure ].We resample linear velocities v    using linear interpolation:v     =   =0       v    , where       are the barycentric weights of the resampled positions x    relative to x    .This preserves the linear velocities of the original element points v   = v    , 0 ≤   <  and preserves the total linear momentum Here we use the linear reproduction property of the resampling interpolation functions N    to equate the barycentric weights of the center of mass with the interpolated barycentric weights of the resample points
Table 3 shows average per-timestep runtime details for several of our examples.For this table, all experiments were run on a desktop with an Intel Core i9-10920X 3.5 GHz LGA 2066 12-Core processor equipped with 32GB RAM.11DISCUSSION AND FUTURE WORKOur hybrid particle/grid iteration improves the state-of-the-art in simulation of volumetric elastic contact.As shown in Section 9.3, our approach improves in conservation and speed over IPC

Table 1 .
Summary of Notation.Number of elements in the one ring (in boundary mesh) of node  x  Mesh node positions at time   v  Mesh node velocities at time   Boundary element nodal linear velocity of element node   m    Resampled boundary element nodal mass of resampled element node   x    Resampled boundary element nodal position of resampled element node   v    Resampled boundary element nodal velocity of resampled element node   Â    Resampled boundary element nodal linear velocity of resampled element node   mi Resampled background grid mass for grid node i xi Resampled element background grid position for grid node i pi Resampled element background grid momentum for grid node i x Barycentric weights of the resampled positions   with respect to original positions   Linear and Angular Momentum Conserving Hybrid Particle/Grid Iteration for Volumetric Elastic Contact Proc.ACM Comput.Graph.Interact.Tech., Vol. 6, No. 3, Article 44.Publication date: August 2023.A