III. Spatial Discretization空间离散化
the discretization of spatially continuous quantities空间的连续离散
a variety of techniques：
- explain the difference between Eulerian and Lagrangian frames of reference欧拉坐标系和拉格朗日坐标系
- grid-based data structures基于网格的数据结构； mesh-free particle-based discretizations基于粒子的无网格离散
- Interpolation techniques
- finite difference 有限差(用于空间导数离散化)
- finite elements有限元(用于空间导数离散化)
3.1 Lagrangian vs. Eulerian
e.g.我们想知道水流速度。方法1：进入水中，看看水通过我们的速度。确定了参考系和测量点，称为固定参考系欧拉参考系（Eulerian reference frames）；方法二:在水中乘坐皮筏，看皮筏的运动速度。参考系随物质运动，称为拉格朗日参考系（Lagrangian reference frames）。 欧拉参考系(固定):常用于流体；使用规则网格和有限差值 fixed reference frames offer many computational advantages 拉格朗日参考系(随材料移动):常用于固体；采用四面体网格或粒子、有限元法。 often used for soft bodies，便于从rest space变形或映射world space计算弹性力 不同的数学标记： 拉，often write the material derivative of a field ??： D y D t \frac{Dy}{Dt} DtDy 应用链用链规则： ? y ? t u ? ▽ y \frac{\partial y}{\partial t} \mathbf u·\bigtriangledown y ?t?y u?▽y，(见2.2.3）
3.2 Grids, Meshes, Particles
Grids
regular grid，常用的空间数据结构，所有边具有相同的长度，即the grid spacing网格间距，常用h或 Δ x \Delta x Δx表示。grid中单独的cube是有8个顶点、12个边、6个面的单元格，可以用一些冗余参数（redundant parameters）来描述:the grid spacing网格间距, the grid resolution分辨率 (i.e. the number of cells in each dimension每一维中有多少格), and the upper and lower extent of the gird. 在空间中是固定的、形状不变的，所以通常采用欧拉坐标系。 在很多实现中，grid是抽象的，每次只有一部分被分配。这种情况下，分层结构（比如octree八叉树或者B-tree）通常被用来代表整个domain（举例：OpenVDB 详见[Museth 2013] ） 仿真变量存储位置？cell centers或者vertices（两种惯例都被广泛应用） 一个特例：staggered grid交错网格（比较少用的叫法是MAC网格，Harlow and Welch [1965]提出），在不同的位置存储不同的变量，常用于流体仿真，使用有限差值时以更小的额外计算量获得二阶精度，代码会比变量co-located的情况复杂。 单个部分的速度 u \mathbf u u作为标量被存储在面的中心，压力被存储在cell centers，如下图。一些作者更喜欢用half-index notation，速度写成 u i − 1 / 2 , j , k \mathbf u_{i-1/2,j,k} ui−1/2,j,k,但这反映在计算机代码中会更不清晰 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-6hFI2O0A-1616637060914)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E4%BA%A4%E9%94%99%E7%BD%91%E6%A0%BC_20200616225803.jpg?v=1&type=image&token=V1:wNUhvpwJLI3OBWzARghYDhJK-TD6Wa_a25701Z7BkMg)]
Meshes
an ill-defined term一个不明确的术语，通常的是指simplicial complexes单纯复形，一个𝑘-simplex包含k+1个all connected相关的顶点，0-simplex就是一个点，1-simplex是线段或边，2-simplex是三角形，3-simplex是四面体，A simplicial complex is decomposes分解 a domain into a set of disjoint不相交的 simplices单形,低维的简单结构组成复杂结构。Automatic tetrahedral meshing自动四面体网格划分仍是难点问题。对任意几何体产生高质量的四面体网格是阻碍它们广泛应用的极大障碍。
Particles
产生tetrahedral meshes四面体网格难度大，转而选择用一系列粒子做几何表示。附近粒子会相互作用，这种方式主要优势是简单，缺点是相互作用会非常复杂（因为空间没有被划分成不相交的集合）
Hybrid Structures 混合结构
combining Lagrangian particles拉格朗日粒子 with regular background grids规则背景网格，比如在the Fluid Implicit Particle (FLIP)流体隐式粒子 和Material Point Method (MPM)物质点法中证明，结合了两种方式的许多优势
3.3 Interpolation
标量场、向量场或张量场的值通常在grid, mesh, or set of particles的离散点上采样，若字段的值不在采样点上，则需要从离散采样点到要求的位置进行插值。
简单的插值举例——单线性插值：
如下图， f ( t ) = ( 1 − t ) f 1 + t f 2 f(t)=(1-t)f_1+tf_2 f(t)=(1−t)f1+tf2 （这里公式可以由 f 2 − f ( t ) 1 − t = f ( t ) − f 1 t \frac{f_2-f(t)}{1-t}=\frac{f(t)-f_1}{t} 1−tf2−f(t)=tf(t)−f1推出，用的较少的一种形式为 f ( t ) = f 1 + t ( f 2 − f 1 ) f(t)=f_1+t(f_2-f_1) f(t)=f1+t(f2−f1)） 在这里我们认为 f ( 0 ) = f 1 , f ( 1 ) = f 2 f(0)=f_1,f(1)=f_2 f(0)=f1,f(1)=f2,计算两终点的加权平均，权值为取样点对面子线段的占比。 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-P1ogHpI8-1616637060915)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E7%AE%80%E5%8D%95%E7%9A%84%E6%8F%92%E5%80%BC%E4%B8%BE%E4%BE%8B_20200617185553.jpg?v=1&type=image&token=V1:ualhrt_RVRmBS70Sksrr-6jjaHBbWGYWPL5OBSbvjd8)] 在任何插值公式中，权值加和为1，形成单位分解（a partition of unity.） 推广generalizing to bilinear, trilinear, and barycentric interpolation
双线性bilinear插值：
如下图，取样点在矩形的四个顶点 f ( s , t ) = ( 1 − s ) ( 1 − t ) f 1 + s ( 1 − t ) f 2 + s t f 3 + ( 1 − s ) t f 4 f(s,t)=(1-s)(1-t)f_1+s(1-t)f_2+stf_3+(1-s)tf_4 f(s,t)=(1−s)(1−t)f1+s(1−t)f2+stf3+(1−s)tf4,权值是取样点相关区域对面区域的占比，可以用两步线性插值推出 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-TLwqRmUC-1616637060917)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E5%8F%8C%E7%BA%BF%E6%80%A7%E6%8F%92%E5%80%BC_20200617193109.jpg?v=1&type=image&token=V1:jujxeatxsJ-vPFrdQRBmzs9GnHNV8vR_RH2y7EGzFaI)]
三线性插值trilinear interpolation：
如下图，取样点在立方体的八个顶点，使用取样点相关体积相对的立方体的占比作为权值，则 f ( s , t , u ) = ( 1 − s ) ( 1 − t ) ( 1 − u ) f 1 + s ( 1 − t ) ( 1 − u ) f 2 + s t ( 1 − u ) f 3 + ( 1 − s ) t ( 1 − u ) f 4 + ( 1 − s ) ( 1 − t ) u f 5 + s ( 1 − t ) u f 6 + s t u f 7 + ( 1 − s ) t u f 8 f(s,t,u)=(1-s)(1-t)(1-u)f_1+s(1-t)(1-u)f_2+st(1-u)f_3+(1-s)t(1-u)f_4+(1-s)(1-t)uf_5+s(1-t)uf_6+stuf_7+(1-s)tuf_8 f(s,t,u)=(1−s)(1−t)(1−u)f1+s(1−t)(1−u)f2+st(1−u)f3+(1−s)t(1−u)f4+(1−s)(1−t)uf5+s(1−t)uf6+stuf7+(1−s)tuf8 可以先双线性插值，使用三个坐标中的两个，然后再在相对的两个面用第三个坐标做单线性插值 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eMwzcG2z-1616637060919)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E4%B8%89%E7%BA%BF%E6%80%A7%E6%8F%92%E5%80%BC_20200617194630.jpg?v=1&type=image&token=V1:9mr6yUkEL2IoET-ek-CUHEr1V5iZceJ0qBb2CjYRIjE)]
重心坐标barycentric coordinates：
在三角形的三个顶点间线性插值， f p = f ( α , β , γ ) = α f a + β f b + γ f c f_p=f(\alpha,\beta,\gamma)=\alpha f_a+\beta f_b+\gamma f_c fp=f(α,β,γ)=αfa+βfb+γfc α , β , γ \alpha,\beta,\gamma α,β,γ是权重，其和为1，p点连接三角形三个顶点，某点权重是对边与p点的子三角形的占比， α = a r e a ( p , b , c ) a r e a ( a , b , c ) \alpha=\frac{area(p,b,c)}{area(a,b,c)} α=area(a,b,c)area(p,b,c), β = a r e a ( p , a , c ) a r e a ( a , b , c ) \beta=\frac{area(p,a,c)}{area(a,b,c)} β=area(a,b,c)area(p,a,c), γ = a r e a ( p , a , b ) a r e a ( a , b , c ) \gamma=\frac{area(p,a,b)}{area(a,b,c)} γ=area(a,b,c)area(p,a,b), 若是在顶点，在该顶点相对的权重为1，其余为0；若在三条边上，则其对面的顶点权重为0，另外两点的权重可看做在该边上的单线性插值 冗余项 γ = 1 − α − β \gamma=1-\alpha-\beta γ=1−α−β, f p = f ( α , β ) = 1 − α − β f a + α f b + β f c = f a + α ( f b − f a ) + β ( f c − f a ) f_p=f(\alpha,\beta)=1-\alpha-\beta f_a+\alpha f_b+\beta f_c=f_a+\alpha(f_b-f_a)+\beta(f_c-f_a) fp=f(α,β)=1−α−βfa+αfb+βfc=fa+