Abstract In recent years, the Smoothed Particle Hydrodynamics (SPH) approach has been increasingly used for multiphase simulations involving interactions between diverse materials. A critical component of an SPH simulator is the pressure solver, which not only facilitates the simulation of compressible or incompressible fluids but also handles contact by preventing penetration between different materials. Currently, most SPH simulations in computer graphics employ implicit dual pressure solvers such as PBF, IISPH, or DFSPH. However, these solvers often exhibit instability when simulating high density ratios. Furthermore, they are difficult to strongly couple with many existing methods for non‐pressure forces, which typically utilize primal formulations. Consequently, pressure and non‐pressure solvers are often only weakly coupled, which can lead to stability issues. We present a novel implicit primal SPH pressure solver designed for multiphase simulations. Our method enables stable simulation of multiple interacting materials with large density ratios. We show that our solver robustly handles ratios of up to 1:1000 (e.g., air‐water interactions) which was not possible with previous implicit SPH pressure solvers. Moreover, we demonstrate how our solver allows for strong coupling with existing implicit simulation methods for viscosity, elasticity, and surface tension. Overall, our strong coupling significantly improves stability in complex multiphase simulations involving fluids, highly viscous materials, and deformable solids.

Abstract Newton's Method is widely used to find the solution of complex non‐linear simulation problems. To guarantee a descent direction, it is common practice to clamp the negative eigenvalues of each element Hessian prior to assembly—a strategy known as Projected Newton (PN)—but this perturbation often hinders convergence. In this work, we observe that projecting only a small subset of element Hessians is sufficient to secure a descent direction. Building on this insight, we introduce Progressively Projected Newton (PPN), a novel variant of Newton's Method that uses the current iterate's residual to cheaply determine the subset of element Hessians to project. The benefit is twofold: most eigendecompositions are avoided and the global Hessian remains closer to its original form, reducing the number of Newton iterations. We compare PPN with PN and Project‐on‐Demand Newton (PDN) in a comprehensive set of experiments covering contact‐free and contact‐rich deformables, co‐dimensional and rigid‐body simulations, and a range of time step sizes, tolerances and resolutions. PPN reduces the amount of element projections in dynamic simulations by one order of magnitude while simultaneously improving convergence, consistently being the fastest solver in our benchmark.

Optimization time integrators are effective at solving complex multi-physics problems including deformable solids with non-linear material models, contact with friction, strain limiting, and so on. For challenging problems, Newton-type optimizers are often used, which necessitates first- and second-order derivatives of the global non-linear objective function. Manually differentiating, implementing, testing, optimizing, and maintaining the resulting code is extremely time-consuming, error-prone, and precludes quick changes to the model, even when using tools that assist with parts of such pipeline. We present SymX, 1 an open source framework that computes the required derivatives of the different energy contributions by symbolic differentiation, generates optimized code, compiles it on-the-fly, and performs the global assembly. The user only has to provide the symbolic expression of each energy for a single representative element in its corresponding discretization and our system will determine the assembled derivatives for the whole simulation. We demonstrate the versatility of SymX in complex simulations featuring different non-linear materials, high-order finite elements, rigid body systems, adaptive discretizations, frictional contact, and coupling of multiple interacting physical systems. SymX’s derivatives offer performance on par with SymPy, an established off-the-shelf symbolic engine, and produces simulations at least one order of magnitude faster than TinyAD, an alternative state-of-the-art integral solution.

Abstract Physics simulation is a cornerstone of many computer graphics applications, ranging from video games and virtual reality to visual effects and computational design. The number of techniques for physically‐based modeling and animation has thus skyrocketed over the past few decades, facilitating the simulation of a wide variety of materials and physical phenomena. This report captures the state‐of‐the‐art of multiphysics simulation for computer graphics applications. Although a lot of work has focused on simulating individual phenomena, here we put an emphasis on methods developed by the computer graphics community for simulating various physical phenomena and materials, as well as the interactions between them. These include combinations of discretization schemes, mathematical modeling frameworks, and coupling techniques. For the most commonly used methods we provide an overview of the state‐of‐the‐art and deliver valuable insights into the various approaches. A selection of software frameworks that offer out‐of‐the‐box multiphysics modeling capabilities is also presented. Finally, we touch on emerging trends in physics‐based animation that affect multiphysics simulation, including machine learning‐based methods which have become increasingly popular in recent years.

Abstract Continuum‐based shell models are an established approach for the simulation of thin deformables in computer graphics. However, existing research in physically‐based animation is mostly focused on shear‐rigid Kirchhoff‐Love shells. In this work we explore three‐director Cosserat (micropolar) shells which introduce additional rotational degrees of freedom. This microrotation field models transverse shearing and in‐plane drilling rotations. We propose an incremental potential formulation of the Cosserat shell dynamics which allows for strong coupling with frictional contact and other physical systems. We evaluate a corresponding finite element discretization for non‐planar shells using second‐order elements which alleviates shear‐locking and permits simulation of curved geometries. Our formulation and the discretization, in particular of the rotational degrees of freedom, is designed to integrate well with typical simulation approaches in physically‐based animation. While the discretization of the rotations requires some care, we demonstrate that they do not pose significant numerical challenges in Newton's method. In our experiments we also show that the codimensional shell model is consistent with the respective three‐dimensional model. We qualitatively compare our formulation with Kirchhoff‐Love shells and demonstrate intriguing use cases for the additional modes of control over dynamic deformations offered by the Cosserat model such as directly prescribing rotations or angular velocities and influencing the shell's curvature.

The numerical simulation of surface tension is an active area of research in many different fields of application and has been attempted using a wide range of methods. Our contribution is the derivation and implementation of an implicit cohesion force based approach for the simulation of surface tension effects using the Smoothed Particle Hydrodynamics (SPH) method. We define a continuous formulation inspired by the properties of surface tension at the molecular scale which is spatially discretized using SPH. An adapted variant of the linearized backward Euler method is used for time discretization, which we also strongly couple with an implicit viscosity model. Finally, we extend our formulation with adhesion forces for interfaces with rigid objects. Existing SPH approaches for surface tension in computer graphics are mostly based on explicit time integration, thereby lacking in stability for challenging settings. We compare our implicit surface tension method to these approaches and further evaluate our model on a wider variety of complex scenarios, showcasing its efficacy and versatility. Among others, these include but are not limited to simulations of a water crown, a dripping faucet, and a droplet toy.

We present a new octree-based neighborhood search method for SPH simulation. A speedup of up to 1.9x is observed in comparison to state-of-the-art methods which rely on uniform grids. While our method focuses on maximizing performance in fixed-radius SPH simulations, we show that it can also be used in scenarios where the particle support radius is not constant thanks to the adaptive nature of the octree acceleration structure. Neighborhood search methods typically consist of an acceleration structure that prunes the space of possible particle neighbor pairs, followed by direct distance comparisons between the remaining particle pairs. Previous works have focused on minimizing the number of comparisons. However, in an effort to minimize the actual computation time, we find that distance comparisons exhibit very high throughput on modern CPUs. By permitting more comparisons than strictly necessary, the time spent on preparing and searching the acceleration structure can be reduced, yielding a net positive speedup. The choice of an octree acceleration structure, instead of the uniform grid typically used in fixed-radius methods, ensures balanced computational tasks. This benefits both parallelism and provides consistently high computational intensity for the distance comparisons. We present a detailed account of high-level considerations that, together with low-level decisions, enable high throughput for performance-critical parts of the algorithm. Finally, we demonstrate the high performance of our algorithm on a number of large-scale fixed-radius SPH benchmarks and show in experiments with a support radius ratio up to 3 that our method is also effective in multi-resolution SPH simulations.

Abstract Visually appealing and vivid simulations of deformable solids represent an important aspect of physically based computer animation. For the temporal discretization, it is customary in computer animation to use first‐order accurate integration methods, such as Backward Euler, due to their simplicity and robustness. Although there is notable research on second‐order methods, their use is not widespread. Many of these well‐known methods have significant drawbacks such as severe numerical damping or scene‐dependent time step restrictions to ensure stability. In this paper, we discuss the most relevant requirements on such methods in computer animation and motivate the interest beyond first‐order accuracy. Keeping these requirements in mind, we investigate several promising methods from the families of diagonally implicit Runge‐Kutta (DIRK) and Rosenbrock methods which currently do not appear to have considerable popularity in this field. We show that the usage of such methods improves the visual quality of physical animations. In addition, we demonstrate that they allow distinctly more control over damping at lower computational cost than classical methods. As part of our theoretical contribution, we review aspects of simulations that are often considered more intricate with higher‐order methods, such as contact handling. To this end, we derive an implicit linearized contact model based on a predictor‐corrector approach that leads to consistent behavior with higher‐order integrators as predictors. Our contact model is well suited for the simulation of stiff, nonlinear materials with the integration methods presented in this paper and more common methods such as Backward Euler alike.

As demands for high-fidelity physics-based animations increase, the need for accurate methods for simulating deformable solids grows. While higherorder finite elements are commonplace in engineering due to their superior approximation properties for many problems, they have gained little traction in the computer graphics community. This may partially be explained by the need for finite element meshes to approximate the highly complex geometry of models used in graphics applications. Due to the additional perelement computational expense of higher-order elements, larger elements are needed, and the error incurred due to the geometry mismatch eradicates the benefits of higher-order discretizations. One solution to this problem is the embedding of the geometry into a coarser finite element mesh. However, to date there is no adequate, practical computational framework that permits the accurate embedding into higher-order elements. We develop a novel, robust quadrature generation method that generates theoretically guaranteed high-quality sub-cell integration rules of arbitrary polynomial accuracy. The number of quadrature points generated is bounded only by the desired degree of the polynomial, independent of the embedded geometry. Additionally, we build on recent work in the Finite Cell Method (FCM) community so as to tackle the severe ill-conditioning caused by partially filled elements by adapting an Additive-Schwarz-based preconditioner so that it is suitable for use with state-of-the-art non-linear material models from the graphics literature. Together these two contributions constitute a general-purpose framework for embedded simulation with higher-order finite elements. We finally demonstrate the benefits of our framework in several scenarios, in which second-order hexahedra and tetrahedra clearly outperform their first-order counterparts.