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.

We present ST-FLIP, a spatiotemporal extension of the Fluid-Implicit Particle (FLIP) method for incompressible free-surface and two-phase liquid simulation. ST-FLIP enables time steps up to an order of magnitude larger than those typically used in CFL-constrained solvers, while preserving detailed flow structures and visual fidelity. It addresses a common failure mode of large time steps in hybrid particle-grid liquid solvers: temporal under-sampling of particle motion produces aliasing-driven free-surface artifacts after projection. Our key idea is to interpret particles as samples in four-dimensional space-time: in addition to standard spatial jittering, we randomize particle positions along the time axis as well and perform particle-to-grid deposition using a separable 4D kernel. This yields a Monte Carlo estimator of perstep time-slab-integrated particle quantities. Although particles are treated as samples in 4D space-time, our approach works as a lightweight plugin by collapsing to slab-integrated 3D grid fields for projection. Building on recent particle-based phase-field work, we reuse the particle-to-grid weight accumulators as a conceptual space-time phase field, providing variable-coefficient projection weights and eliminating the need for per-step surface reconstruction. The method can be easily integrated into existing FLIP/PIC or APIC solvers with negligible additional computational cost per time step. The effectiveness of our approach is demonstrated through a series of comparisons with state-of-the-art solvers, yielding several-fold speedups for multi-billion-particle simulations at high effective 3D resolutions on a single workstation.

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.

Neural shape representation generally refers to representing 3D geometry using neural networks, e.g., computing a signed distance or occupancy value at a specific spatial position. In this paper we present a neural-network architecture suitable for accurate encoding of 3D shapes in a single forward pass. Our architecture is based on a multi-scale hybrid system incorporating graph-based and voxel-based components, as well as a continuously differentiable decoder. The hybrid system includes a novel way of voxelizing point-based features in neural networks by projecting the point "feature-field" onto a grid. This projection is insensitive to local point density, and we show that it can be used to obtain smoother and more detailed reconstructions, in particular when combined with oriented point clouds as input. Our architecture also requires only a single forward pass, instead of the latent-code optimization used in auto-decoder methods. Furthermore, our network is trained to solve the well-established eikonal equation and only requires knowledge of the zero-level set for training and inference. We additionally propose a modification to the aforementioned loss function for the case that surface normals are not well defined, e.g., in the context of non-watertight surfaces and non-manifold geometry. Overall, our method consistently outperforms other baselines on the surface reconstruction task across a wide variety of datasets, while being more computationally efficient and requiring fewer parameters.

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.

We present a novel implicit porous flow solver using SPH, which maintains fluid incompressibility and is able to model a wide range of scenarios, driven by strongly coupled solid-fluid interaction forces. Many previous SPH porous flow methods reduce particle volumes as they transition across the solid-fluid interface, resulting in significant stability issues. We instead allow fluid and solid to overlap by deriving a new density estimation. This further allows us to extend SPH pressure solvers to take local porosity into account and results in strict enforcement of incompressibility. As a result, we can simulate porous flow using physically consistent pressure forces between fluid and solid. In contrast to previous SPH porous flow methods, which use explicit forces for internal fluid flow, we employ implicit non-pressure forces. These we solve as a linear system and strongly couple with fluid viscosity and solid elasticity. We capture the most common effects observed in porous flow, namely drag, buoyancy and capillary action due to adhesion. To achieve elastic behavior change based on local fluid saturation, such as bloating or softening, we propose an extension to the elasticity model. We demonstrate the efficacy of our model with various simulations that showcase the different aspects of porous flow behavior. To summarize, our system of strongly coupled non-pressure forces and enforced incompressibility across overlapping phases allows us to naturally model and stably simulate complex porous interactions.

Capturing the visually compelling features of large-scale water phenomena, such as the spray clouds of crashing waves, stormy seas, or waterfalls, involves simulating not only the water but also the motion of the air interacting with it. However, current solutions in the visual effects industry still largely rely on single-phase solvers and non-physical "white-water" heuristics. To address these limitations, we present Phase-Field-FLIP (PF-FLIP), a hybrid Eulerian/Lagrangian method for the fully physics-based simulation of very large-scale, highly turbulent multiphase flows at high Reynolds numbers and high fluid density contrasts. PF-FLIP transports mass and momentum in a consistent, non-dissipative manner and, unlike most existing multiphase approaches, does not require a surface reconstruction step. Furthermore, we employ spatial adaptivity across all critical components of the simulation algorithm, including the pressure Poisson solver. We augment PF-FLIP with a dual multiresolution scheme that couples an efficient treeless adaptive grid with adaptive particles, along with a fast adaptive Poisson solver tailored for high-density-contrast multiphase flows. Our method enables the simulation of two-phase flow scenarios with a level of physical realism and detail previously unattainable in graphics, supporting billions of particles and adaptive 3D resolutions with thousands of grid cells per dimension on a single workstation.

The use of simulation in robotics is increasingly widespread for the purpose of testing, synthetic data generation and skill learning. A relevant aspect of simulation for a variety of robot applications is physics-based simulation of robot-object interactions. This involves the challenge of accurately modeling and implementing different mechanical systems such as rigid and deformable bodies as well as their interactions via constraints, contact or friction. Most state-of-the-art physics engines commonly used in robotics either cannot couple deformable and rigid bodies in the same framework, lack important systems such as cloth or shells, have stability issues in complex friction-dominated setups or cannot robustly prevent penetrations. In this paper, we propose a framework for strongly coupled simulation of rigid and deformable bodies with focus on usability, stability, robustness and easy access to state-of-the-art deformation and frictional contact models. Our system uses the Finite Element Method (FEM) to model deformable solids, the Incremental Potential Contact (IPC) approach for frictional contact and a robust second order optimizer to ensure stable and penetration-free solutions to tight tolerances. It is a general purpose framework, not tied to a particular use case such as grasping or learning, it is written in C++ and comes with a Python interface. We demonstrate our system’s ability to reproduce complex real-world experiments where a mobile vacuum robot interacts with a towel on different floor types and towel geometries. Our system is able to reproduce 100% of the qualitative outcomes observed in the laboratory environment. The simulation pipeline, named Stark (the German word for strong, as in strong coupling) is made open-source.

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.

Abstract We present a strongly coupled method for the robust simulation of linear magnetic rigid bodies. Our approach describes the magnetic effects as part of an incremental potential function. This potential is inserted into the reformulation of the equations of motion for rigid bodies as an optimization problem. For handling collision and friction, we lean on the Incremental Potential Contact (IPC) method. Furthermore, we provide a novel, hybrid explicit / implicit time integration scheme for the magnetic potential based on a distance criterion. This reduces the fill‐in of the energy Hessian in cases where the change in magnetic potential energy is small, leading to a simulation speedup without compromising the stability of the system. The resulting system yields a strongly coupled method for the robust simulation of magnetic effects. We showcase the robustness in theory by analyzing the behavior of the magnetic attraction against the contact resolution. Furthermore, we display stability in practice by simulating exceedingly strong and arbitrarily shaped magnets. The results are free of artifacts like bouncing for time step sizes larger than with the equivalent weakly coupled approach. Finally, we showcase the utility of our method in different scenarios with complex joints and numerous magnets.

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.

Dynamics simulation with frictional contacts is important for a wide range of applications, from cloth simulation to object manipulation. Recent methods using smoothed lagged friction forces have enabled robust and differentiable simulation of elastodynamics with friction. However, the resulting frictional behavior can be inaccurate and may not converge to analytic solutions. Here we evaluate the accuracy of lagged friction models in comparison with implicit frictional contact systems. We show that major inaccuracies near the stick-slip threshold in such systems are caused by lagging of friction forces rather than by smoothing the Coulomb friction curve. Furthermore, we demonstrate how systems involving implicit or lagged friction can be correctly used with higher-order time integration and highlight limitations in earlier attempts. We demonstrate how to exploit forward-mode automatic differentiation to simplify and, in some cases, improve the performance of the inexact Newton method. Finally, we show that other complex phenomena can also be simulated effectively while maintaining smoothness of the entire system. We extend our method to exhibit stick-slip frictional behavior and preserve volume on compressible and nearly-incompressible media using soft constraints.

Physical simulation is an indispensable component of robotics simulation platforms that serves as the basis for a plethora of research directions. Looking strictly at robotics, the common characteristic of the most popular physics engines, such as ODE, DART, MuJoCo, bullet, SimBody, PhysX or RaiSim, is that they focus on the solution of articulated rigid bodies with collisions and contacts problems, while paying less attention to other physical phenomena. This restriction limits the range of addressable simulation problems, rendering applications such as soft robotics, cloth simulation, simulation of viscoelastic materials, and fluid dynamics, especially surface swimming, infeasible. In this work, we present Gazebo Fluids, an open-source extension of the popular Gazebo robotics simulator that enables the interaction of articulated rigid body dynamics with particle-based fluid and deformable solid simulation. We implement fluid dynamics and highly viscous and elastic material simulation capabilities based on the Smoothed Particle Hydrodynamics method. We demonstrate the practical impact of this extension for previously infeasible application scenarios in a series of experiments, showcasing one of the first self-propelled robot swimming simulations with SPH in a robotics simulator.

Abstract Throughout the past decades, the graphics community has spent major resources on the research and development of physics simulators on the mission to computer‐generate behaviors achieving outstanding visual effects or to make the virtual world indistinguishable from reality. The variety and impact of recent research based on Smoothed Particle Hydrodynamics (SPH) demonstrates the concept's importance as one of the most versatile tools for the simulation of fluids and solids. With this survey, we offer an overview of the developments and still‐active research on physics simulation methodologies based on SPH that has not been addressed in previous SPH surveys. Following an introduction about typical SPH discretization techniques, we provide an overview over the most used incompressibility solvers and present novel insights regarding their relation and conditional equivalence. The survey further covers recent advances in implicit and particle‐based boundary handling and sampling techniques. While SPH is best known in the context of fluid simulation we discuss modern concepts to augment the range of simulatable physical characteristics including turbulence, highly viscous matter, deformable solids, as well as rigid body contact handling. Besides the purely numerical approaches, simulation techniques aided by machine learning are on the rise. Thus, the survey discusses recent data‐driven approaches and the impact of differentiable solvers on artist control. Finally, we provide context for discussion by outlining existing problems and opportunities to open up new research directions.

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.

Iterative solvers are widely used to accurately simulate physical systems. These solvers require initial guesses to generate a sequence of improving approximate solutions. In this contribution, we introduce a novel method to accelerate iterative solvers for rod dynamics with graph networks (GNs) by predicting the initial guesses to reduce the number of iterations. Unlike existing methods that aim to learn physical systems in an end-to-end manner, our approach guarantees long-term stability and therefore leads to more accurate solutions. Furthermore, our method improves the run time performance of traditional iterative solvers for rod dynamics. To explore our method we make use of position-based dynamics (PBD) as a common solver for physical systems and evaluate it by simulating the dynamics of elastic rods. Our approach is able to generalize across different initial conditions, discretizations, and realistic material properties. We demonstrate that it also performs well when taking discontinuous effects into account such as collisions between individual rods. Finally, to illustrate the scalability of our approach, we simulate complex 3D tree models composed of over a thousand individual branch segments swaying in wind fields.

We propose a novel implicit density projection approach for hybrid Eulerian/Lagrangian methods like FLIP and APIC to enforce volume conservation of incompressible liquids. Our approach is able to robustly recover from highly degenerate configurations and incorporates volume-conserving boundary handling. A problem of the standard divergence-free pressure solver is that it only has a differential view on density changes. Numerical volume errors, which occur due to large time steps and the limited accuracy of pressure projections, are invisible to the solver and cannot be corrected. Moreover, these errors accumulate over time and can lead to drastic volume changes, especially in long-running simulations or interactive scenarios. Therefore, we introduce a novel method that enforces constant density throughout the fluid. The density itself is tracked via the particles of the hybrid Eulerian/Lagrangian simulation algorithm. To achieve constant density, we use the continuous mass conservation law to derive a pressure Poisson equation which also takes density deviations into account. It can be discretized with standard approaches and easily implemented into existing code by extending the regular pressure solver. Our method enables us to relax the strict time step and solver accuracy requirements of a regular solver, leading to significantly higher performance. Moreover, our approach is able to push fluid particles out of solid obstacles without losing volume and generates more uniform particle distributions, which makes frequent particle resampling unnecessary. We compare the proposed method to standard FLIP and APIC and to previous volume correction approaches in several simulations and demonstrate significant improvements in terms of incompressibility, visual realism, and computational performance.

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.

In this article, we present a novel method for the robust handling of static and dynamic rigid boundaries in Smoothed Particle Hydrodynamics (SPH) simulations. We build upon the ideas of the density maps approach which has been introduced recently by Koschier and Bender. They precompute the density contributions of solid boundaries and store them on a spatial grid which can be efficiently queried during runtime. This alleviates the problems of commonly used boundary particles, like bumpy surfaces and inaccurate pressure forces near boundaries. Our method is based on a similar concept but we precompute the volume contribution of the boundary geometry. This maintains all benefits of density maps but offers a variety of advantages which are demonstrated in several experiments. First, in contrast to the density maps method we can compute derivatives in the standard SPH manner by differentiating the kernel function. This results in smooth pressure forces, even for lower map resolutions, such that precomputation times and memory requirements are reduced by more than two orders of magnitude compared to density maps. Furthermore, this directly fits into the SPH concept so that volume maps can be seamlessly combined with existing SPH methods. Finally, the kernel function is not baked into the map such that the same volume map can be used with different kernels. This is especially useful when we want to incorporate common surface tension or viscosity methods that use different kernels than the fluid simulation.

In this paper, we present a novel method for the robust handling of static and dynamic rigid boundaries in Smoothed Particle Hydrodynamics (SPH) simulations. We build upon the ideas of the density maps approach which has been introduced recently by Koschier and Bender. They precompute the density contributions of solid boundaries and store them on a spatial grid which can be efficiently queried during runtime. This alleviates the problems of commonly used boundary particles, like bumpy surfaces and inaccurate pressure forces near boundaries. Our method is based on a similar concept but we precompute the volume contribution of the boundary geometry and store it on a grid. This maintains all benefits of density maps but offers a variety of advantages which are demonstrated in several experiments. Firstly, in contrast to the density maps method we can compute derivatives in the standard SPH manner by differentiating the kernel function. This results in smooth pressure forces, even for lower map resolutions, such that precomputation times and memory requirements are reduced by more than two orders of magnitude compared to density maps. Furthermore, this directly fits into the SPH concept so that volume maps can be seamlessly combined with existing SPH methods. Finally, the kernel function is not baked into the map such that the same volume map can be used with different kernels. This is especially useful when we want to incorporate common surface tension or viscosity methods that use different kernels than the fluid simulation.

We present a strong fluid-rigid coupling for Smoothed Particle Hydrodynamics (SPH) fluids and rigid bodies with particle-sampled surfaces. The approach interlinks the iterative pressure update at fluid particles with a second SPH solver that computes artificial pressure at rigid-body particles. The introduced SPH rigid-body solver models rigid-rigid contacts as artificial density deviations at rigid-body particles. The corresponding pressure is iteratively computed by solving a global formulation that is particularly useful for large numbers of rigid-rigid contacts. Compared to previous SPH coupling methods, the proposed concept stabilizes the fluid-rigid interface handling. It significantly reduces the computation times of SPH fluid simulations by enabling larger time steps. Performance gain factors of up to 58 compared to previous methods are presented. We illustrate the flexibility of the presented fluid-rigid coupling by integrating it into DFSPH, IISPH, and a recent SPH solver for highly viscous fluids. We further show its applicability to a recent SPH solver for elastic objects. Large scenarios with up to 90 M particles of various interacting materials and complex contact geometries with up to 90 k rigid-rigid contacts are shown. We demonstrate the competitiveness of our proposed rigid-body solver by comparing it to Bullet.

In this paper we introduce a novel micropolar material model for the simulation of turbulent inviscid fluids. The governing equations are solved by using the concept of Smoothed Particle Hydrodynamics (SPH). As already investigated in previous works, SPH fluid simulations suffer from numerical diffusion which leads to a lower vorticity, a loss in turbulent details and finally in less realistic results. To solve this problem we propose a micropolar fluid model. The micropolar fluid model is a generalization of the classical Navier-Stokes equations, which are typically used in computer graphics to simulate fluids. In contrast to the classical Navier-Stokes model, micropolar fluids have a microstructure and therefore consider the rotational motion of fluid particles. In addition to the linear velocity field these fluids also have a field of microrotation which represents existing vortices and provides a source for new ones. However, classical micropolar materials are viscous and the translational and the rotational motion are coupled in a dissipative way. Since our goal is to simulate turbulent fluids, we introduce a novel modified micropolar material for inviscid fluids with a non-dissipative coupling. Our model can generate realistic turbulences, is linear and angular momentum conserving, can be easily integrated in existing SPH simulation methods and its computational overhead is negligible. Another important visual feature of turbulent liquids is foam. Therefore, we present a post-processing method which considers microrotation in the foam particle generation. It works completely automatic and requires only one user-defined parameter to control the amount of foam.

Abstract In this paper, we present a novel direct solver for the efficient simulation of stiff, inextensible elastic rods within the position‐based dynamics (PBD) framework. It is based on the XPBD algorithm, which extends PBD to simulate elastic objects with physically meaningful material parameters. XPBD approximates an implicit Euler integration and solves the system of non‐linear equations using a non‐linear Gauss–Seidel solver. However, this solver requires many iterations to converge for complex models and if convergence is not reached, the material becomes too soft. In contrast, we use Newton iterations in combination with our direct solver to solve the non‐linear equations which significantly improves convergence by solving all constraints of an acyclic structure (tree), simultaneously. Our solver only requires a few Newton iterations to achieve high stiffness and inextensibility. We model inextensible rods and trees using rigid segments connected by constraints. Bending and twisting constraints are derived from the well‐established Cosserat model. The high performance of our solver is demonstrated in highly realistic simulations of rods consisting of multiple 10 000 segments. In summary, our method allows the efficient simulation of stiff rods in the PBD framework with a speedup of two orders of magnitude compared to the original XPBD approach.

Abstract In this paper we present a novel operator splitting approach for corotated FEM simulations. The deformation energy of the corotated linear material model consists of two additive terms. The first term models stretching in the individual spatial directions and the second term describes resistance to volume changes. By formulating the backward Euler time integration scheme as an optimization problem, we show that the first term is invariant to rotations. This allows us to use an operator splitting approach and to solve both terms individually with different numerical methods. The stretching part is solved accurately with an optimization integrator, which can be done very efficiently because the system matrix is constant over time such that its Cholesky factorization can be precomputed. The volume term is solved approximately by using the compliant constraints method and Gauss‐Seidel iterations. Further, we introduce the analytic polar decomposition which allows us to speed up the extraction of the rotational part of the deformation gradient and to recover inverted elements. Finally, this results in an extremely fast and robust simulation method with high visual quality that outperforms standard corotated FEMs by more than two orders of magnitude and even the fast but inaccurate PBD and shape matching methods by more than one order of magnitude without having their typical drawbacks. This enables a very efficient simulation of complex scenes containing more than a million elements.

Abstract In this paper, we present a novel physically consistent implicit solver for the simulation of highly viscous fluids using the Smoothed Particle Hydrodynamics (SPH) formalism. Our method is the result of a theoretical and practical in‐depth analysis of the most recent implicit SPH solvers for viscous materials. Based on our findings, we developed a list of requirements that are vital to produce a realistic motion of a viscous fluid. These essential requirements include momentum conservation, a physically meaningful behavior under temporal and spatial refinement, the absence of ghost forces induced by spurious viscosities and the ability to reproduce complex physical effects that can be observed in nature. On the basis of several theoretical analyses, quantitative academic comparisons and complex visual experiments we show that none of the recent approaches is able to satisfy all requirements. In contrast, our proposed method meets all demands and therefore produces realistic animations in highly complex scenarios. We demonstrate that our solver outperforms former approaches in terms of physical accuracy and memory consumption while it is comparable in terms of computational performance. In addition to the implicit viscosity solver, we present a method to simulate melting objects. Therefore, we generalize the viscosity model to a spatially varying viscosity field and provide an SPH discretization of the heat equation.

In this paper we present a robust remeshing-free cutting algorithm on the basis of the eXtended Finite Element Method (XFEM) and fully implicit time integration. One of the most crucial points of the XFEM is that integrals over discontinuous polynomials have to be computed on subdomains of the polyhedral elements. Most existing approaches construct a cut-aligned auxiliary mesh for integration. In contrast, we propose a cutting algorithm that includes the construction of specialized quadrature rules for each dissected element without the requirement to explicitly represent the arising subdomains. Moreover, we solve the problem of ill-conditioned or even numerically singular solver matrices during time integration using a novel algorithm that constrains non-contributing degrees of freedom (DOFs) and introduce a preconditioner that efficiently reuses the constructed quadrature weights. Our method is particularly suitable for fine structural cutting as it decouples the added number of DOFs from the cut's geometry and correctly preserves geometry and physical properties by accurate integration. Due to the implicit time integration these fine features can still be simulated robustly using large time steps. As opposed to this, the vast majority of existing approaches either use remeshing or element duplication. Remeshing based methods are able to correctly preserve physical quantities but strongly couple cut geometry and mesh resolution leading to an unnecessary large number of additional DOFs. Element duplication based approaches keep the number of additional DOFs small but fail at correct conservation of mass and stiffness properties. We verify consistency and robustness of our approach on simple and reproducible academic examples while stability and applicability are demonstrated in large scenarios with complex and fine structural cutting.

In this paper we present a novel Smoothed Particle Hydrodynamics (SPH) method for the efficient and stable simulation of incompressible fluids. The most efficient SPH-based approaches enforce incompressibility either on position or velocity level. However, the continuity equation for incompressible flow demands to maintain a constant density and a divergence-free velocity field. We propose a combination of two novel implicit pressure solvers enforcing both a low volume compression as well as a divergence-free velocity field. While a compression-free fluid is essential for realistic physical behavior, a divergence-free velocity field drastically reduces the number of required solver iterations and increases the stability of the simulation significantly. Thanks to the improved stability, our method can handle larger time steps than previous approaches. This results in a substantial performance gain since the computationally expensive neighborhood search has to be performed less frequently. Moreover, we introduce a third optional implicit solver to simulate highly viscous fluids which seamlessly integrates into our solver framework. Our implicit viscosity solver produces realistic results while introducing almost no numerical damping. We demonstrate the efficiency, robustness and scalability of our method in a variety of complex simulations including scenarios with millions of turbulent particles or highly viscous materials.

In this paper we present an hp-adaptive algorithm to generate discrete higher-order polynomial Signed Distance Fields (SDFs) on axis-aligned hexahedral grids from manifold polygonal input meshes. Using an orthonormal polynomial basis, we efficiently fit the polynomials to the underlying signed distance function on each cell. The proposed error-driven construction algorithm is globally adaptive and iteratively refines the SDFs using either spatial subdivision ( h-refinement) following an octree scheme or by cell-wise adaption of the polynomial approximation's degree ( p-refinement). We further introduce a novel decision criterion based on an error-estimator in order to decide whether to apply p- or h-refinement. We demonstrate that our method is able to construct more accurate SDFs at significantly lower memory consumption compared to previous approaches. While the cell-wise polynomial approximation will result in highly accurate SDFs, it can not be guaranteed that the piecewise approximation is continuous over cell interfaces. Therefore, we propose an optimization-based post-processing step in order to weakly enforce continuity. Finally, we apply our generated SDFs as collision detector to the physically-based simulation of geometrically highly complex solid objects in order to demonstrate the practical relevance and applicability of our method.

We present a novel algorithm to extract the rotational part of an arbitrary 3 X 3 matrix. This problem lies at the core of two popular simulation methods in computer graphics, the co-rotational Finite Element Method and Shape Matching techniques. In contrast to the traditional method based on polar decomposition, degenerate configurations and inversions are handled robustly and do not have to be treated in a special way. In addition, our method can be implemented with only a few lines of code without branches which makes it particularly well suited for GPU-based applications. We demonstrate the robustness, coherence and efficiency of our method by comparing it to stabilized polar decomposition in several simulation scenarios.

We present a new method for particle based fluid simulation, using a combination of Projective Dynamics and Smoothed Particle Hydrodynamics (SPH). The Projective Dynamics framework allows the fast simulation of a wide range of constraints. It offers great stability through its implicit time integration scheme and is parallelizable in large parts, so that it can make use of modern multi core CPUs. Yet existing work only uses Projective Dynamics to simulate various kinds of soft bodies and cloth. We are the first ones to incorporate fluid simulation into the Projective Dynamics framework. Our proposed fluid constraints are derived from SPH and seamlessly integrate into the existing method. Furthermore, we adapt the solver to handle the constantly changing constraints that appear in fluid simulation. We employ a highly parallel matrix-free conjugate gradient solver, and thus do not require expensive matrix factorizations.

Abstract We propose a position‐based approach for large‐scale simulations of rigid bodies at interactive frame rates. Our method solves positional constraints between rigid bodies and can therefore be seamlessly integrated into other position‐based methods. Interaction of particles and rigid bodies through common constraints enables two‐way coupling with deformables. The method exhibits exceptional performance and stability while being user controllable and easy to implement. Various results demonstrate the practicability of our method for the resolution of collisions, contacts, stacking and joint constraints. Copyright © 2014 John Wiley & Sons, Ltd.

Abstract Volumetric scalar data sets are common in many scientific, engineering and medical applications where they originate from measurements or simulations. Furthermore, they can represent geometric scene content, e.g. as distance or density fields. Often isosurfaces are extracted, either for indirect volume visualization in the former category, or to simply obtain a polygonal representation in case of the latter. However, even moderately sized volume data sets can result in complex isosurfaces which are challenging to recompute in real time, e.g. when the user modifies the isovalue or when the data itself are dynamic. In this paper, we present a GPU‐friendly algorithm for the extraction of isosurfaces, which provides adaptive level of detail rendering with view‐dependent tessellation. It is based on a longest edge bisection scheme where the resulting tetrahedral cells are subdivided into four hexahedra, which then form the domain for the subsequent isosurface extraction step. Our algorithm generates meshes with good triangle quality even for highly non‐linear scalar data. In contrast to previous methods, it does not require any stitching between regions of different levels of detail. As all computation is performed at run time and no pre‐processing is required, the algorithm naturally supports dynamic data and allows us to change isovalues at any time.

Abstract Interactive rigid body simulation is an important part of many modern computer tools, which no authoring tool nor game engine can do without. Such high‐performance computer tools open up new possibilities for changing how designers, engineers, modelers and animators work with their design problems. This paper is a self contained state‐of‐the‐art report on the physics, the models, the numerical methods and the algorithms used in interactive rigid body simulation all of which have evolved and matured over the past 20 years. Furthermore, the paper communicates the mathematical and theoretical details in a pedagogical manner. This paper is not only a stake in the sand on what has been done, it also seeks to give the reader deeper insights to help guide their future research .

Abstract The dynamic simulation of mechanical effects has a long history in computer graphics. The classical methods in this field discretize Newton's second law in a variety of Lagrangian or Eulerian ways, and formulate forces appropriate for each mechanical effect: joints for rigid bodies; stretching, shearing or bending for deformable bodies and pressure, or viscosity for fluids, to mention just a few. In the last years, the class of position‐based methods has become popular in the graphics community. These kinds of methods are fast, stable and controllable which make them well‐suited for use in interactive environments. Position‐based methods are not as accurate as force‐based methods in general but they provide visual plausibility. Therefore, the main application areas of these approaches are virtual reality, computer games and special effects in movies. This state‐of‐the‐art report covers the large variety of position‐based methods that were developed in the field of physically based simulation. We will introduce the concept of position‐based dynamics, present dynamic simulation based on shape matching and discuss data‐driven upsampling approaches. Furthermore, we will present several applications for these methods.

Computing ambient occlusion in screen-space (SSAO) is a common technique in real-time rendering applications which use rasterization to process 3D triangle data. However, one of the most critical problems emerging in screen-space is the lack of information regarding occluded geometry which does not pass the depth test and is therefore not resident in the G-buffer. These occluded fragments may have an impact on the proximity-based shadowing outcome of the ambient occlusion pass. This not only decreases image quality but also prevents the application of SSAO on multiple layers of transparent surfaces where the shadow contribution depends on opacity. We propose a novel approach to the SSAO concept by taking advantage of per-pixel fragment lists to store multiple geometric layers of the scene in the G-buffer, thus allowing order independent transparency (OIT) in combination with high quality, opacity-based ambient occlusion (OITAO). This A-buffer concept is also used to enhance overall ambient occlusion quality by providing stable results for low-frequency details in dynamic scenes. Furthermore, a flexible compression-based optimization strategy is introduced to improve performance while maintaining high quality results.

Abstract We present graphics processing unit (GPU) data structures and algorithms to efficiently solve sparse linear systems that are typically required in simulations of multi‐body systems and deformable bodies. Thereby, we introduce an efficient sparse matrix data structure that can handle arbitrary sparsity patterns and outperforms current state‐of‐the‐art implementations for sparse matrix vector multiplication. Moreover, an efficient method to construct global matrices on the GPU is presented where hundreds of thousands of individual element contributions are assembled in a few milliseconds. A finite‐element‐based method for the simulation of deformable solids as well as an impulse‐based method for rigid bodies are introduced in order to demonstrate the advantages of the novel data structures and algorithms. These applications share the characteristic that a major computational effort consists of building and solving systems of linear equations in every time step. Our solving method results in a speed‐up factor of up to 13 in comparison to other GPU methods.

Abstract This paper describes an impulse‐based dynamic simulation method for articulated bodies which has a linear time complexity. Existing linear‐time methods are either based on a reduced‐coordinate formulation or on Lagrange multipliers. The impulse‐based simulation has advantages over these well‐known methods. Unlike reduced‐coordinate methods, it handles nonholonomic constraints like velocity‐dependent ones and is very easy to implement. In contrast to Lagrange multiplier methods the impulse‐based approach has no drift problem and an additional stabilisation is not necessary. The presented method computes a simulation step inO(n) time for acyclic multi‐body systems containing equality constraints. Closed kinematic chains can be handled by dividing the model into different acyclic parts. Each of these parts is solved independently from each other. The dependencies between the single parts are solved by an iterative method. In the same way inequality constraints can be integrated in the simulation process in order to handle collisions and permanent contacts with dynamic and static friction. Copyright © 2007 John Wiley & Sons, Ltd.