Fluid Simulation

100 particles colliding with sphere

Building off of our cloth simulator, we set up a particle-based fluid simulator. Because cloth is simulated by a system of point masses and springs, by removing the springs, and rendering a sphere for each point mass, we adapt the system to model liquid. Each particle has the properties mass, pressure, density, position, forces, and velocity. We implemented a 3D smoothed particle hydrodynamics (SPH) solver to calculate the positions of the point masses as a result of the accelerations due to density, pressure, and the external forces (gravity) acting on it. At a high level, our simulation is governed by the following physical equations referenced in Penn State EGEE 520:

$$\frac{dv_i}{dt} = a_i^{pressure} + a_i^{viscosity} + a_i^{gravity}$$ $$a_i^{pressure} = \langle -\frac{1}{\text{density}}\nabla \text{pressure}\rangle_i$$ $$a_i^{viscosity} = \langle -\frac{\text{viscosity}}{\text{density}}\nabla \nabla \text{velocity}\rangle_i$$ $$a_i^{gravity}\approx [0\ 0\ -g]^{T}$$

These equations are known as the Navier-Stokes equations, which model the motion of viscous fluids. Since these equations are continuous (and therefore unable to be computed with 100% accuracy), we approximate them using summations - more specifically, given particle $p_i$, we iterate over $p_i$'s neighbors, (calculated using the algorithm outlined below), weighting each neighbor $p_j$ according to a kernel function $K(r)$. This kernel assigns higher values to neighbors that are closer spatially to $p_i$, and vice versa. Summing all neighboring contributions with respect to pressure and viscosity (pressure relying on the gradient $\nabla K(r)$, and viscosity relying on the second derivative $\frac{d^2}{dr^2}K(r)$) gives us $p_i$'s respective accelerations. Finally, with this particle's acceleration and mass, we can calculate the forces acting on $p_i$ for a single timestep $\Delta t$.

In order to create a arc-like stream of fluid originating from the teapot spout and ending at the teacup, we initialize each particle with a starting velocity. It took many trials and errors to fine tune the $x$, $y$, and $z$ components of our starting velocity to achieve the perfect movement. To simulate a continuous stream, we repeatedly spawn batches of particles at the spout every ~100 time steps.

STL Renderer

Tea set made in blender

Rendered tea set STL

Using Blender, we modeled a teacup and imported Blender's premade teapot mesh to create our scene. We positioned the teacup and teapot next to each other with the teapot slightly raised above it and rotated, as if it were pouring tea into the cup. We exported this as an STL file which stores the values of the points and normals of the triangles in the triangle mesh. Given this information, we built an STL renderer to render the teaset model.

Collisions

Particles colliding with the teacup with normal shading

Liquid pooling in the teacup with mirror shading

Liquid pooling in the teacup with Phong shading

Having loaded the STL file, we store the triangles as a bounding volume hierarchy (BVH) tree to check for intersections between the particles and triangles efficiently. Thus, in addition to spheres, planes, and boxes, we support collisions with the triangle meshes that make up the teacup and teapot.

Optimizations

To iterate over each particle's neighbors for the SPH solver, testing proximity for each particle was computationally infeasible. Thus, the acceleration structure we chose to tackle this problem was a grid-based hash map, modeled after Fast and Efficient Nearest Neighbor Search for Particle Simulations. When computing the effect of nearby particles on density, pressure, and acceleration, we only consider particles in the same grid cell—approximating the kernel defined for SPH.

Similarly, our BVH tree optimizes our system to handle collisions by only checking for particle-triangle collisions in a bounding box, based on the position of the particle. Because we recursively check for collisions with a bounding box, the tree structure limits the depth of recursion.

Unlike light rays, which can intersect both the left & right bounding boxes, a single point mass can only occupy at most one of the two at any given time. This allows us to short circuit our traversal sooner and optimize the runtime of intersection tests.

We also utilized OpenMP to parallelize certain tasks. To do this, we #pragma omp parallel for directive above each iterative loop, with an additional keyword of collapse(n) for n-level nested loops (for example, n=3 with marching cubes algorithms).

Surface Reconstruction

To accomplish surface reconstruction, we used the Marching Cubes algorithm. At a high level, we takes the following steps:

1. Initialize a grid of vertices across our world space (representing the vertices of our Marching Cubes)
2. For each point mass, identify the 8 nearest vertices and add itself to their lists of neighbors
3. Compute the density at each vertex
4. Compute the *activation* at each vertex (1 if $\rho>\rho_{\text{rest}}$, 0 otherwise)
5. For each cube, use the cube's 8 vertex *activations* to access the corresponding triangulated cube from a table of all 256 triangulated cubes
6. Render the triangles of each selected triangulated cube onto the scene

Triangulated Cubes

Shaders

Tea set with mirror shading, plane with floral texture, and translucent water

In general, for rendering, an object uses the most recently-bound shader. Specifically, for HW4's "Texture" shader and plane collision object, after applying the texture to the cloth (whose $uv$ coordinates were calculated correctly), it is implicitly applied to the plane. Since, by default, the plane doesn't set any $uv$ attribute, the texture is incorrectly mapped. By modifying plane.cpp such that the bounds of the image span the plane, the texture is correctly mapped. For our purposes, this results in a glossy marble table surface for our teaset to rest on!

Mirror shading with dining room cube map

Closeup of mirror shading with dining room cube map

Honey-like tea with sky cube map

---