WENO Reconstruction Background In many numerical tasks, a continuous function is represented by discrete samples and reconstructed when values between sample locations are needed. In this post, we consider a function $f(x)$ defined on the interval $[0, L]$. We divide the interval into $N$ uniform subintervals with grid spacing $$ \Delta x = \frac{L}{N}. $$ The grid points are located at $$ x_i = i\Delta x, \qquad i = 0, 1, \ldots, N, $$

Problem Description In this blog, we want to talk about an algorithm which can compute prefix sum (or pre-sum) of an array in parallel. Formally, for an input array a[N], we want an output array c[N] satisfying : $$ c[i] = \sum_{0\leq j < i} a[j] $$ Since it doesn’t contain the value a[i], we call it “exclusive pre-sum”. Fenwick Tree In sequential implementation, this problem is realy easy and intuitive.

Some Other Things for Implicit MPM In the last blog, we talked about some basic things for implicit MPM. But there are also some details we should mention in this continued blog. Hessian Matrix In the last blog, we defined $E(\mathbf{u}) = \frac{1}{2}(\mathbf{u}-\mathbf{u^*})\mathbf{M}(\mathbf{u}-\mathbf{u^*}) + \Phi(\mathbf{x}^n+\Delta t \mathbf{u})$ as the energy for a system, our goal is to minimize it. We know the gradient of it is : $$ \nabla_\mathbf{u}E = \mathbf{M}(\mathbf{u}-\mathbf{u}^*) + \Delta t \frac{\partial \Phi}{\partial\mathbf{x}}|_{(\mathbf{x}^n+\Delta t \mathbf{u})} = \mathbf{M}(\mathbf{u}-\mathbf{u}^*) - \Delta t \mathbf{f}(\mathbf{x}^n+\Delta t \mathbf{u}) $$

Implicit Update Scheme for MPM The main purpose of this article is to talk about how to apply implicit update in MPM framework. To know more details about MPM, you can read This SIGGRAPH 2016 Course Notes and A Material Point Method For Snow Simulation (Alexey Stomakhin et al., 2013) . Framework of MPM Currently, I prefer to use APIC scheme to transfer velocity between particles and grids. The framework I always use is listed here:

Conjugate Gradient Method in Solving Linear Equation System Linear Equation System In fluid simulation, sometimes we need to solve a Linear Equation System like doing projection in Eulerian fluid. A linear Equation System can be written as a matrix form: $$ Ax = b $$ where $A$ is a $N\times N$ matrix and $x,b$ are $N\times 1$ vectors. Actually, we have a lot of methods to solve this problem, like Gauss-Seidel iteration, Gaussian elimination, etc.

Finite Elements Method for Simulating Elastic Body What do we do in simulating spring-mass system? Recall what we did in simulating Spring-Mass system: Firstly, we describe deformation of springs using $\mathbf{F} = \frac{\mathbf{x}_1 - \mathbf{x}_2}{l_0}$. Something noticable here is that $\mathbf{F}$ is a vector and it is related to rotation of springs. Then, we define strain $G = |\mathbf{F}| - 1$. The purpose of $G$ is that it is a scale and only related to stretch and contraction of springs.