论文笔记_SIGGRAPH2019会前课程：An Introduction to Physics-Based Animation_2

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.我们想知道水流速度。方法一：进入水中看水经过我们的速度。其参考系和测量点是确定的，这样的固定参考系叫欧拉参考系（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 𝑦：$\frac{Dy}{Dt}$
欧，应用链式法则：$\frac{\partial y}{\partial t}+\mathbf{u}\cdot ▽y$，（将固定点的变化率和通过流场的平流分开，见2.2.3）

3.2 Grids, Meshes, Particles

Grids

regular grid，常用的空间数据结构，所有边具有相同的长度，即the grid spacing网格间距，常用h或$\Delta 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的情况复杂。
单个部分的速度$\mathbf{u}$作为标量被存储在面的中心，压力被存储在cell centers，如下图。一些作者更喜欢用half-index notation，速度写成${\mathbf{u}}_{i-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$
（这里公式可以由$\frac{f_2-f(t)}{1-t}=\frac{f(t)-f_1}{t}$推出，用的较少的一种形式为$f(t)=f_1+t(f_2-f_1)$）
在这里我们认为$f(0)=f_1,f(1)=f_2$,计算两终点的加权平均，权值为取样点对面子线段的占比。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(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+st{f}_3+(1-s)t{f}_4$,权值是取样点相关区域对面区域的占比，可以用两步线性插值推出
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(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+st(1-u)f_3+(1-s)t(1-u)f_4+(1-s)(1-t)u{f}_5+s(1-t)u{f}_6+stu{f}_7+(1-s)tu{f}_8$
可以先双线性插值，使用三个坐标中的两个，然后再在相对的两个面用第三个坐标做单线性插值
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(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(\alpha ,\beta ,\gamma )=\alpha f_a+\beta f_b+\gamma f_c$
$\alpha ,\beta ,\gamma$是权重，其和为1，p点连接三角形三个顶点，某点权重是对边与p点的子三角形的占比，
$\alpha =\frac{area(p,b,c)}{area(a,b,c)}$,$\beta =\frac{area(p,a,c)}{area(a,b,c)}$,$\gamma =\frac{area(p,a,b)}{area(a,b,c)}$,
若是在顶点，在该顶点相对的权重为1，其余为0；若在三条边上，则其对面的顶点权重为0，另外两点的权重可看做在该边上的单线性插值
冗余项$\gamma =1-\alpha -\beta$,
$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)$
这种方式往往会产生更高效的代码
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-877NRwVr-1616637060919)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E9%87%8D%E5%BF%83%E5%9D%90%E6%A0%87_20200617201714.jpg?v=1&type=image&token=V1:i_iviXBL2AOWLbehKS3P71QKXbPKxXgpNuS4aHvYvOA)]
上面的区域是带符号的区域，约定逆时针方向为正，顺时针方向为负。
重心坐标可以为负但其和为1，也可用来判断点是否在三角形内（全为正则在三角形内，有一到两个为负则在三角形外）
寻找重心坐标的另一方法见[Marschner and Shirley 2015]

Polynomial interpolation.多项式插值

𝑛 + 1 distinct data points$(t_n,f_n)$，polynomial多项式 of degree 𝑛,拉格朗日插值：
$f(t)={\sum }_{i=1}^{i=1}{f}_i{l}_i(t)$
$l_i(t)=\frac{{\prod }_{k\ne i}(t-t_k)}{{\prod }_{k\ne i}(t_i-t_k)}$
多项式也可以用monomial basis elements定义，导出the ill-conditioned Vandermonde matrix范德蒙矩阵，或者使用牛顿插值。详见[Heath 2002]（textbook）

Approximating functions函数逼近

插值函数经过一组离散样本可以定义样本间的函数值，但近似函数不是必须经过数据点，近似一组样本数据的最常用方法：least squares最小二乘 interpolation
data is noisy or a system is overdetermined,点插值结果可能不理想，找一个含未知参数的符合数据的模型，最小化模型与数据间的误差。
最小二乘中，给定m个数据点$(x_i,f_i)$,寻找数据模型，一组n个基函数${\varphi }_i(x)$，n个系数${\alpha }_i$
得到近似函数$\varphi (x)={\alpha }_1{\varphi }_1(x)+...+{\alpha }_n{\varphi }_n(x)$
方差(使其最小)：${\sum }_{i=1}^m|\varphi (x_i)-f_i{|}^2$
基函数选定后使系数拟合数据
最优化：$argmi{n}_{\alpha }||A\alpha -\mathbf{f}||$
$A=\left[\begin{array}{ccc}{\varphi }_1(x_1)& {\varphi }_2(x_1)......& {\varphi }_n(x_1)\\ {\varphi }_1(x_2)& {\varphi }_2(x_2)......& {\varphi }_n(x_2)\\ ...& & \\ {\varphi }_1(x_m)& {\varphi }_2(x_m)......& {\varphi }_n(x_m)\end{array}\right]$
$\alpha =\left[\begin{array}{c}{\alpha }_1\\ ...\\ {\alpha }_n\end{array}\right]$,$\mathbf{f}=\left[\begin{array}{c}{f}_1\\ ...\\ f_n\end{array}\right]$
正规方程（normal equations）:$A^{T}A\alpha =A^T\mathbf{f}$
可用于解非奇异最小二乘问题（poorly conditioned，有更好替代见[Heath 2002].）
带权方差${\sum }_{i=1}^m{\omega }_i|\varphi (x_i)-f_i{|}^2$
权值可以表示每个误差项在总误差中的贡献程度
推广：这些权重可以在定义域上变化，比如moving least squares approaches移动最小二乘法可用于离散数据，从非结构化数据插值重构函数，见[Shen et al.2005]
piecewise polynomial curves分段多项式曲线（比如贝赛尔曲线、B-splines B样条）的近似有时也被用来根据样本数据重构函数，见[Marschner and Shirley 2015]
人物模型的几何变形，有专用于插值的坐标系，能使用更普遍的结构。比如mean-value coordinates中值坐标（见[Floater 2003; Ju et al. 2005]）；harmonic coordinate（见[Joshi et al. 2007]）；green coordinates（见[Lipman et al. 2008]）；bounded biharmonic weights（见[Jacobson et al. 2011]）

3.4 Finite Differences

are commonly used to discretize the differential
operators for numerical solution.常用于离散数值解的微分算子

前向差分

回忆微分的定义，无穷小量换成很小但有限的数h，h
越小近似的误差越小
$\frac{d}{dx}f(x)\approx \frac{f(x+h)-f(x)}{h}$
h和误差的依赖关系可以使用函数f在点x处的泰勒级数（在x附近收敛）来研究，特别地，
$f(x+h)=f(x)+h{f}^{\prime }(x)+...+\frac{h^n}{n!}{f}^{(n)}(x)+...$
除以h，并整理，得
$f^{\prime }(x)=\frac{f(x+h)-f(x)}{h}-\frac{h}{2}{f}^{\prime \prime }(x)+\frac{h^2}{6}{f}^{\prime \prime \prime }(x)+...+\frac{h^{n-1}}{n!}{f}^{(n)}(x)+...$
即$f^{\prime }(x)=\frac{f(x+h)-f(x)}{h}+O(h)$
这种有限差分对$f^{\prime }(x)$的近似叫前向差分，主要误差项为O(h)，是一阶精确的。

后向差分：

$f^{\prime }(x)\approx \frac{f(x)-f(x-h)}{h}$
函数f在点x处的泰勒级数展开：
$f(x-h)=f(x)-h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{\prime \prime }(x)-...+\frac{h^n}{n!}{f}^{(n)}(x)+...$
整理得 $f^{\prime }(x)=\frac{f(x)-f(x-h)}{h}+\frac{h}{2}{f}^{\prime \prime }(x)-\frac{h^2}{6}{f}^{\prime \prime \prime }(x)+...+\frac{h^{n-1}}{n!}{f}^{(n)}(x)+...$
即$f^{\prime }(x)=\frac{f(x)-f(x-h)}{h}+O(h)$
后向差分也是一阶精确
二者主要误差项的区别是符号。一阶精确度的意义是，如果量取样点中间的间隔h是二分的，主要误差项也会被二分，as it is linear in the size of the mesh spacing

central difference中心差分

常见的二阶精度有限差分
$f^{\prime }(x)\approx \frac{f(x+h)-f(x-h)}{2h}$
泰勒展开
$f(x+h)=f(x)+h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{\prime \prime }(x)+\frac{h^3}{6}{f}^{\prime \prime \prime }(x)+O(h^4)$
$f(x-h)=f(x)-h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{\prime \prime }(x)-\frac{h^3}{6}{f}^{\prime \prime \prime }(x)+O(h^4)$
两式做差并除以2h(只剩偶数项)
$\frac{f(x+h)-f(x-h)}{2h}=f^{\prime }(x)+\frac{h^2}{6}{f}^{\prime \prime \prime }(x)+O(h^3)=f^{\prime }(x)+O(h^2)$
主要误差项是O(h^2),二阶精确。如果把网格间距从h降到h/2，误差项会缩1/4，因此，近似地，随着h降低，误差以一种更高阶的方式降低更快。
在选择有限差分策略时，精度的阶数只是几个考虑因素之一。其他因素有：overall stability of the scheme，the nature of the errors produced by the scheme(e.g., dissipative or dispersive),the conservation properties of the scheme(比如当离散偏微分方程中的平流项advective terms时，中心差分法结合certain time discretization methods会导致不好的震荡或不稳定性)
Upwind discretizations逆风离散 常用于advection(热的水平对流；空气的水平运动)，使用local流向来选择前向或后向差分，或更高阶的单侧离散。
用有限差分近似逼近高阶导数 e.g.用中心差分近似二阶导
$f^{\prime \prime }(x)=\frac{f^{\prime }(x+\frac{h}{2})-f^{\prime }(x-\frac{h}{2})}{h}+O(h^2)$
$=\frac{\frac{f(x+h)-f(x)}{h}-\frac{f(x)-f(x-h)}{h}}{h}+O(h^2)$
$=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}+O(h^2)$
泰勒展开分析，the odd errors terms of 𝑓 (𝑥 − ℎ) and 𝑓 (𝑥 + ℎ) 完全抵消，主要误差项分子上剩下$O(h^4)$，所以它是二阶精度的。

Laplacian operator.

考虑去近似 函数u的拉普拉斯算子，即使用有限差分方法近似$\Delta u$
二维的情况，$\Delta u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}$,这两个二阶导都可以用中心差分的策略近似(在上边的例子中就是)
$\frac{{\partial }^{2}u(x,y)}{\partial x^2}\approx \frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2}$
$\frac{{\partial }^{2}u(x,y)}{\partial y^2}\approx \frac{u(x,y+h)-2u(x,y)+u(x,y-h)}{h^2}$
两式相加
$\frac{{\partial }^{2}u(x,y)}{\partial x^2}+\frac{{\partial }^{2}u(x,y)}{\partial y^2}$
$\approx \frac{u(x+h,y)+u(x-h,y)+u(x,y+h)+u(x,y-h)-4u(x,y)}{h^2}$
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ksew7ryz-1616637060921)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E4%BA%8C%E7%BB%B4%E7%BD%91%E6%A0%BC_20200621021916.jpg?v=1&type=image&token=V1:5eVHG1HwO5DcegWlBfEVV3Fq24keKptDceA6allfkHQ)]
如图，周围的取样点被用于有限差分近似，is called the stencil漏字板，模板 of the finite difference scheme
$\Delta u(x,y)\approx \frac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}}{h^2}$
这个Laplacian所绘的特殊模板非常常用，被称为五点模板（five-point stencil），类似的，一维——三点，二维——七点。
离散一维的泊松方程：$u_{xx}=f,x\in \Omega$
在大小为8的一维网格上,在网格每个点进行有限差分（$f^{\prime \prime }(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2}+O(h^2)$）
将Laplacian在网格上离散，得到一个稀疏但全局耦合所有变量的矩阵。因为微分运算符提供函数的local information，其有限差分离散化会得到稀疏矩阵
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wQ57y85w-1616637060922)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E4%B8%80%E7%BB%B4grid_20200621023934.jpg?v=1&type=image&token=V1:Di1-Jt_BmkR_02Q2FvFHjRuWVQdpCuFoCaMNOrier9g)]
$\frac{u_{i-1}-2u_i+u_{i+1}}{h^2}=f_i$
将八个点处的该方程组合成一个线性系统。两端点仅有一边邻居，不满足stencil的要求。数学上，该方程的解需要确定边界条件，边界条件用于在边界点上离散方程。我们假设使用狄利克雷边界条件$u(x)=\bar{u}(x),x\in \partial \Omega$
则$\frac{1}{h^2}\left[\begin{array}{cccccc}-2& 1& 0& 0& 0& 0\\ 1& -2& 1& 0& 0& 0\\ 0& 1& -2& 1& 0& 0\\ 0& 0& 0& 1& -2& 1\\ 0& 0& 0& 0& 1& -2\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right]=\left[\begin{array}{c}{f}_1-\frac{\bar{u}(x_0)}{h^2}\\ f_2\\ f_3\\ f_4\\ f_5\\ f_6-\frac{\bar{u}(x_7)}{h^2}\end{array}\right]$
每行对应三个内部样本，第一行和最后一行只有两个非零项，第三项是边界项而非未知数，因此出现在方程右边。
Neumann诺依曼 boundary conditions指定导数都存在。若设边界导数值为0，$\frac{\partial u(x)}{\partial x}=0\in \partial \Omega$
就像我们处理固体边界上的压强一样，离散的拉普拉斯方程就变成了线性方程组
$\frac{1}{h^2}\left[\begin{array}{cccccc}-1& 1& 0& 0& 0& 0\\ 1& -2& 1& 0& 0& 0\\ 0& 1& -2& 1& 0& 0\\ 0& 0& 0& 1& -2& 1\\ 0& 0& 0& 0& 1& -1\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right]=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\\ f_5\\ f_6\end{array}\right]$
第一行和最后一行只有两个非零项，对角线项减少了，对第一行和最后一行，第三项等于解，矩阵项消掉了1

3.5 Finite Elements

严谨的数学推导有些复杂，需要引入任意测试函数和分部积分。但是广泛用于图形学的其基本思想和线性有限元其实很简单。
主要思想：离散函数空间，可以在这个有限的空间里表示和解决问题。将一个object分成元素的有限集，元素不相交且提供object的空间离散化.公共元素是简单元素（2D-三角形，3D-四面体）、变形的hypercubes（2D-四边形，3D-六面体）。然后在这些元素上定义基函数，最常见的基础函数空间是分段线性的，只讨论这种情况。
待解的方程中的任意函数被投影到这个分段线性空间（Galerkin projection），给定空间离散化，我们找到与感兴趣的函数最接近的分段线性函数。在这个有限空间中解这个问题，得到一个分段线性解，是原始问题的近似解。

Application to Soft Bodies.

有限元常用于elastic bodies.我们使用一个分段线性基函数来表示变形函数$\mathbf{x}(\mathbf{u})$,这个函数的梯度是分段常数,所以要计算变形梯度$\mathbf{F}$，在一个element上是常量。
in the deformed element，计算中心坐标
$\mathbf{x}={\mathbf{x}}_0+\alpha ({\mathbf{x}}_1-{\mathbf{x}}_0)+\beta ({\mathbf{x}}_2-{\mathbf{x}}_0)$
令$u_{ab}$为向量$u_a-u_b$,类似$\mathbf{x}$,写成矩阵形式：
$\mathbf{u}={\mathbf{u}}_0+({\mathbf{u}}_{10}$ ${\mathbf{u}}_{20})$$(\alpha$ $\beta {)}^T$
$\mathbf{x}={\mathbf{x}}_0+({\mathbf{x}}_{10}$ ${\mathbf{x}}_{20})$$(\alpha$ $\beta {)}^T$
矩阵的列向量是沿着元素边缘的向量,
用$\mathbf{u}$表示$(\alpha$ $\beta {)}^T$带入x的方程，得deformation function
$\mathbf{x}(\mathbf{u})={\mathbf{x}}_0+({\mathbf{x}}_{10}$ ${\mathbf{x}}_{20})({\mathbf{u}}_{10}$ ${\mathbf{u}}_{20}{)}^{-1}(\mathbf{u}-{\mathbf{u}}_0)$
函数的梯度关于$\mathbf{u}$,变形梯度：
$\mathbf{F}=\frac{\partial \mathbf{x}}{\partial \mathbf{u}}=({\mathbf{x}}_{10}$ ${\mathbf{x}}_{20})({\mathbf{u}}_{10}$ ${\mathbf{u}}_{20}{)}^{-1}$
通过组合向量沿着元素的边组成的矩阵，使我们可计算单个元素的变形梯度，对一个矩阵求逆然后进行矩阵乘法运算。有了变形梯度，我们就可以计算任何我们想要的应力值，由此计算弹性力。
若对象的静止形状不变，待求逆的矩阵是常量，可提前计算。如下图。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-IUFwf6Jw-1616637060924)(3.%E7%A9%BA%E9%97%B4%E7%A6%BB%E6%95%A3%E5%8C%96_md_files/%E9%9D%99%E6%AD%A2%E5%BD%A2%E7%8A%B6%E4%B8%8D%E5%8F%98_20200621221048.jpg?v=1&type=image&token=V1:vjXnt8AeN0JBi8sc5BIGRsvvYCyqQiVPJeS_rEQ2zpQ)]
在2.2.2中提到的strain model应变模型co-rotated model（可能是最常见的）中，计算极分解$\mathbf{F}=\mathbf{Q}F$,
计算strain:$ϵ=\frac{1}{2}(F+F^T)-\mathbf{I}$
计算stress:$\sigma =\lambda Tr(ϵ)\mathbf{I}+2\mu ϵ$
力的得出：$\mathbf{f}=\mathbf{Q}\sigma {\mathbf{n}}_i$
$n_i$是in the rest configuration,力作用点对面的边(2D)或面(3D)的区域加权法线
算法伪码（与质点弹簧的伪码稍有不同）

 1：for Particle p : particles do
 2: 	p.frc = 0
 3: 	p.frc += p.mass*gravity
 4: end for
 5: for Element e : elements do
 6: 	Matrix3x3 F = Matrix3x3 (x1-x0, x2-x0, x3-x0)* inverse(Matrix3x3(u1-u0, u2-u0, u3-u0))
 7: 	PolarDecomp (F, Q, Ftilde)
 8: 	Matrix3x3 strain = 1/2 * (Ftilde + transpose (Ftilde)) - I
 9: 	Matrix3x3 stress = lambda * trace(strain) * I + 2 * mu * strain
10: 	for 𝑖 = 0 to 3 do
11: 		particles[e.node[i]].frc += Q * stress * e.normal[i]
12:		end for
13: end for
14: for Particle p : particles do
15: 	p.vel += dt*(p.frc / p.mass)
16: 	p.pos += dt*(p.vel)
17: end for

V. Temporal Discretization时间离散化

• the explicit form of temporal discretization andintegration时间离散与积分的显式形式
  • 同样的 truncation error截断误差，两种方法会有不同的表现：the trapezoidal rule for first-order ordinary differential equations vs. the midpoint method for first-order ordinary differential equations一阶常微分方程的梯形法则 和 一阶常微分方程的中点法
  • 二阶常微分方程的显示积分-symplectic Euler method辛欧拉方法（lower order of accuracy精度）
• implicit integration隐式积分（couples together a system of partial differential equations耦合偏微分方程组）
  the linearly implicit Euler method线性显式欧拉方法：by explicitlylinearizing the equations of motion通过将运动方程显式线性化
  briefly discuss alternatives for solving the resulting linear system and the underlyingnonlinear system简要讨论求解得到的线性方程组和潜在的非线性方程组的可替代方法

仿真效果随时间变化，时间离散化是为了以固定频率输出画面。
timestep大小可变，选择固定的大小以简化、避免不期望的材质表现的变化。
一般地，画面的每一帧都有几个timesteps，但现在的仿真系统致力于减少每帧timesteps

4.1 Explicit Integration

时间$t+\Delta t$更新的状态仅由时间t计算的量来描述：
${\mathbf{x}}_p(t+\Delta t)={\mathbf{x}}_p(t)+\Delta t\cdot \mathbf{v}({\mathbf{x}}_p,t)$
这是一个显示积分器。
不同的性能有不同的积分技巧，看起来非常相似的积分技巧可能在实践中有很大不同的表现。

4.1.1 Trapezoidal Rule梯形法则 vs. Midpoint Method中点法

梯形法则

类似于积分求曲边梯形面积，计算质点位置处的速度，假设以该速度移动了质点整个的timestep，计算末速度，二者均值为计算质点位移的速度：
$\mathbf{x}(t+\Delta t)=\mathbf{x}(t)+\frac{\Delta t}{2}(\mathbf{v}(\mathbf{x}(t),t)+\mathbf{v}(\mathbf{x}+\Delta t\mathbf{v}(\mathbf{x},t),t)$

中点法

计算质点位置处的速度,假设以该速度移动了质点一半的timestep，计算该速度，后者为计算质点位移的速度：
$\mathbf{x}(t+\Delta t)=\mathbf{x}(t)+\Delta t(\mathbf{v}(\mathbf{x}+\frac{\Delta t}{2}\mathbf{v}(\mathbf{x},t),t))$

都是second order Runge-Kutta (RK2) methods；
（使用泰勒级数分析）二阶精度；local truncation errors截断误差 of $O(\Delta t^3)$
都常用（在上下文中使用欧拉积分提高精度值得额外计算时）
二者表现大不相同：梯形法则求均值，更平滑，但有时会出现过度阻尼；中点法更容易受噪声或错认信号的影响
数值方法比它的准确性更重要，在分析过程中看起来很微妙的差异在实践中可能很显著

4.1.2 Symplectic Euler（Section 1有介绍）

对于二阶微分方程，比如牛顿第二定律，三种明显的数值方法来更新粒子位置
${\mathbf{x}}_p(t+\Delta t)={\mathbf{x}}_p(t)+dt\cdot {\mathbf{v}}_p(t)$
${\mathbf{x}}_p(t+\Delta t)={\mathbf{x}}_p(t)+\frac{dt}{2}\cdot ({\mathbf{v}}_p(t)+{\mathbf{v}}_p(t+\Delta t))$
——sometimes called Improved Euler（second-order accurate）
${\mathbf{x}}_p(t+\Delta t)={\mathbf{x}}_p(t)+dt\cdot {\mathbf{v}}_p(t+\Delta t)$
——Symplectic Euler
特别地，for pure elasticity problems，improved Euler
is unconditionally unstable for open-ended (i.e. infinite) time integration无限时间差值，即无论timestep有多小都会有分歧，实践中更常用辛欧拉(尽管是低阶精度)，再次证明数值方法的选择不仅仅是简单的精度顺序，而且通常需要数值实验来确定最佳选择。

4.2 Implicit Integration

对于Stiff problem刚性问题（材料有很强的抗变形能力），显示积分器要求timestep尽可能小来保持稳定，对于零长度弹簧的简单例子，$\mathbf{f}=-k\mathbf{x}$，初位置${\mathbf{x}}_0$,初速度0，质量m，使用辛欧拉积分器经过一个timestep：$\mathbf{x}(\Delta t)=(1-\frac{\Delta t^{2}k}{m}){\mathbf{x}}_0$,如果$\Delta t>\sqrt{\frac{m}{k}}$会越过原点,如果$\Delta t>\sqrt{\frac{2m}{k}}$弹簧会伸长，解会发散。k越大timestep就要越小以保持稳定，为了避免timestep的局限，转向隐式积分。
显式辛欧拉积分器：
$\mathbf{v}(t+\Delta t)=\mathbf{v}(t)+\Delta t\cdot {\mathbf{M}}^{-1}\mathbf{f}(\mathbf{x}(t),t)$
$\mathbf{x}(t+\Delta t)=\mathbf{x}(t)+\Delta t\cdot \mathbf{v}(t+\Delta t)$
隐式形式：
$\mathbf{v}(t+\Delta t)=\mathbf{v}(t)+\Delta t\cdot {\mathbf{M}}^{-1}\mathbf{f}(\mathbf{x}(t+\Delta t),t+\Delta t)$
$\mathbf{x}(t+\Delta t)=\mathbf{x}(t)+\Delta t\cdot \mathbf{v}(t+\Delta t)$
区别：隐式积分中力使用$t+\Delta t$计算，新的状态无法直接算出，但符合该等式。这是复杂的非线性问题，牛顿方法经常与各种性能增强方法一起使用。
考虑这个问题最简单的情况，我们将系统围绕当前状态线性化，并解决线性化问题（a single Newton iteration）
对soft bodies采用非线性动力系统：
$\mathbf{K}(\mathbf{x}-{\mathbf{x}}_0)+\mathbf{D}(\dot{x})+\mathbf{M}\ddot{x}={\mathbf{f}}_{ext}$
${\mathbf{x}}_0$是静息状态，将其线性化得到：
$\mathbf{K}\mathbf{x}-\mathbf{K}{\mathbf{x}}_0+\mathbf{D}(\dot{x})+\mathbf{M}\ddot{x}={\mathbf{f}}_{ext}$
合力：
${\mathbf{f}}_{ext}-\mathbf{K}\mathbf{x}+\mathbf{K}{\mathbf{x}}_0-\mathbf{D}(\dot{x})$
将这些力及隐式积分的第二个式子代入隐式积分的第一个式子：
$\mathbf{v}(t+\Delta t)=\mathbf{v}(t)+\Delta t\cdot {\mathbf{M}}^{-1}[{\mathbf{f}}_{ext}-\mathbf{K}(\mathbf{x}(t)+\Delta t\cdot \mathbf{v}(t+\Delta t))+\mathbf{K}{\mathbf{x}}_0-\mathbf{D}\mathbf{v}(t+\Delta t)]$
整理以得到线性系统
$(\mathbf{M}+\Delta t^2\mathbf{K}+\Delta t\mathbf{D})\mathbf{v}(t+\Delta t)=\mathbf{M}\mathbf{v}(t)+\Delta t(-\mathbf{K}(\mathbf{x}(t)-{\mathbf{x}}_0)+{\mathbf{f}}_{ext})$
——线性隐式欧拉，可从一阶泰勒展开推导出来，右边对应一个显式的忽略阻尼力的欧拉步骤，左边的矩阵在这个显式的欧拉步骤中充当一个平滑滤波器。利用左边矩阵的稀疏结构，采用预条件共轭梯度对该系统进行求解。该积分器也是用牛顿法求解非线性问题的完全隐式向后欧拉积分器的第一步，因此也被称为semi-implicit Euler半隐式欧拉或者semi-implicit Backward Euler半隐式后向欧拉。