GAMES 图形学系列笔记（三十七）

GAMES201：高级物理引擎实战指南2020 - P4：Lecture 4 欧拉视角 - GAMES-Webinar - BV1ZK411H7Hc

那么欢迎大家再次回到我们的高级物理引擎，实战课程，我们已经到第四讲了，不知不觉这个课程已经过去差不多一半了，那么今天呢我们讲一讲欧拉视角下的流体模拟。

ok那么在开始讲之前，还是忍不住分享一下论坛上面的一些，大家提交的新的作业，感觉还是很有意思啊，有同学写了这个三维的fm，用上了新的viewer啊，有同学写了这个分型嗯。

有同学呢写了这个vertex ring lip morgan，应该还是一个非常有技术含量的作业，ok那么对这些就是homework林的一些新的作品，我觉得大家可以去论坛上面把代码下下来，自己跑一跑。

然后学习学习呃，这些东西是怎么写的，觉得还是很有意思的，看看别人代码的时候能理解很多新的东西。

好然后homework肯定呢我们当时发了一个写的好讲，然后最近中科院物理所还写了一篇科普文章，大概是讲卡门涡街的这个标题叫做醉酒的蝴蝶，飞不出花花世界的原因，找到了，我也不太知道是什么意思。

但是反正是个讲流体的文章，然后里面用到的demo，其实就是我们和默克林有一个同学的作业，大家可以去看看，这个文章写得还是挺有意思的，主要是讲各种雷洛树下的这个残留啊，卡门涡街还有残留的行为。

非常通俗易懂，然后有很多同学问我啊，homework一之前要求是实现影视，实现这个时间计分器和显示时间积分器的对比，但是有些sober没法实现显示时间的计算器，咋办呢，啊这个其实不用担心。

你可以只只实现影视时间计算器，然后在影视时间计分器内部进行对比，比如说你可以实现不同的影视时间，积分器的格式，比如说jacob啊，你可以或者写这个country gradients啊。

或者mt grade都可以，这些新的格式呢我们今天会介绍到，ok然后论坛呢我们最近来了一次功能，升级到太极，升级到了v零点06：12，然后论坛有了一次很大的更新，可以呃支持贴公式，我试了一下。

还是非常非常实用的功能，然后也可以贴一些呃，贴些图，然后效果还是非常好的，非常感谢，work class，同学帮我们进行了这个非常辛苦的升级工作，现在论坛应该可用性啊比之前要好了很多，然后我们第五讲。

我们呃专门讲多体问题与我方法，我们请了著名的张星星师兄来给我们讲这方面，因为他是这方面的专家，那一讲呢我们会提到多体问题，和他们与柏松问题的联系呃，然后还会张欣怡同学还会给大家分享，分享我方法的乐趣。

然后从直观的角度，引导同学们认识几种不同的快速求和方法呃。

为了欢迎我们的张鑫鑫师兄。

我还特地自己也写了一个窝方法的sober，但是基于他的之前的一个c加加程序写的，我这边可以给大家演示一下，是一个vertex ring lip morgan的程序，ok，还是很有意思。

这个代码只有几十行，大家可以哦，对这个这个代码我还没分享到论坛上，大家可以下载下来看一下，当然也是基本上可以实时的跑的。

ok看这个代码还是，就90行捕捞的代码，然后非常容易理解，下一讲啊，就专门请了我们的客座讲师张艺兴师兄，给我们讲一讲呃窝方法。

然后呢7月6号我们打算空一周，这样给大家多留一些时间去实现一下，开放作业一呃，7月11号呢我们开放作业一评选截止，然后我们会做一些点评，然后发出一些嗯选出一些很好的作品，然后也会发出一些奖品啊。

这个有同学反应好像比较卡，那么呃没事，这个课程这一讲结束以后，我会把程序上传到论坛上面，大家可以自己在自己的机器上面跑，那肯定就会肯定就不卡了，那么课程纪念品呢，我们homework 0的获奖同学。

我们已经全部都寄到了，然后呃当然请大家查收，因为有可能有些地方会比较慢一些，有些地方比较快一些，有些地方可能还没有击倒，然后大家请留意自己的呃短信啊，邮箱啊，然后尽可能的呃确保不要被。

不要被别人拿走了之类的，这个非常这边非常感谢负责后勤的东西，帮我们定制了这么多呃杯子，然后负责邮寄工作，这个真的是一个非常辛苦的工作，对这个有同学问这是不是杯子，它当然是一个杯子呃。

不是后面的这个纪念品，不是后面的这个电源充电器啊，是这个杯子，还有杯盖，这个应该比较容易看出来呃，ok今天内容还是内容还是挺多的，我就不多说废话了，之前我们主要讲的是拉格朗日视角。

那么记得拉格朗日视角是什么，拉格朗日视角这些例子是随波逐流的，他们会随这里的材料在游戏里面移动，那么拉格朗日视角的例子，往往会用来记录一下，到底自己的位置和自己的速度在哪，或者会记录一些温度啊。

或者其他的属性对吧，那么它的关键就是我们的sensor，或者这个传感器呃，你可以想象成这个小的流体的粒子，它是随着你的流体或者固体在空间中移动的，今天呢我们讲一讲这个欧拉视角。

那么欧拉视角呢大家记住它的特点是什么，它的特点是你的这个采样的点，在空间中是完全不动的，ok有同学说我这个非常卡，我是不是应该把我这边的对我，因为我开了一个直播的同时，我嗯我在网上又开了一个网页。

开了一个虎牙的网页来看我自己的直播，可能这样会导致比较卡，我现在把我自己的那个网页关了，但这样的问题就是说我可能会看不到呃，弹幕或者提问，那么如果有问题的话，大家可以在课程群里面发发一微信群。

这样我不知道这样会不会不卡一些，试一试，今天嗯回到这个欧拉视角，欧拉视角它的特点是，它的采样点是永远不不会移动的，然后它的采样点回答什么问题，回答问题是在我这个点这边。

我的呃材料它以什么样的速度在穿过我自己，ok那么这个是欧拉视角的一个直观的理解，你可以认为就是很多在你在水里面插了很多桩，这些桩永远不动，但是这个妆每一个妆啊，都是知道自己周围的水流流速的好。

那么这个是图片部分了，我们下面进入稍微复杂一些的公式部分，呃今天呢我们会稍微讲一点点数学，但是不会讲特别多，然后这门课呢主要还是侧重于直观上的推导，呃我们偶尔会提到一些有限体积，有限微分啊。

有限语言之类的基础知识，但是我们只会提到非常最为直观的一部分，我们不会讲特别深的这些lumerical scheme的理论，如果大家想深入学习欧拉流体模拟啊，可以看一看这本书，这本书写的非常好。

我自己也是看这本书入门的，叫做free simulation for computer graphics，它的作者是robert bryson，这个是一个大佬级的人物了。

在流体计算机通信学里面的流体模拟里面呃，基本上是无人不知无人不晓的一个人物，但是顺便说一句，张星星师兄就是这位大佬的学生，名师出高徒呀，然后呃今天这讲呢，其实我们只是快速的过一下这个。

欧拉法流体模拟的一些内容，ok那么我们先overview一下，看一看这个欧拉法流体模拟到底是在干啥，首先我们要介绍一个概念，这个概念叫做material relatives，或者叫做材料导数。

中文的材料导数这个翻译应该没有什么问题，嗯他其实呢联系了欧拉视角和拉格尔视角，这个台词老师是什么意思呢，它的符号一般是用大d大写的d，然后而不是我们的偏微分符号或者小写的d。

然后呃材料导数它其实有两个部分，当你去求一个物理量的大d by d t的时候，他其实呃它有两个成分，第一个成分是这个物理量关于时间的导数，就partial by partial t对吧。

第二个成分呢是呃，材料在这点速度点乘上一个grad，然后直观上的什么，这个人是什么意思呢，直观上来就说你这个物理量，由于呃你是呃这个是material问题，那也就是说这个采样点。

是随着你的材料一起在移动的，那你这个彩随着材料一起移动，采样点，它的物理上的变化有两个成分对吧，第一个成分呢是你假设你即使不动，那你的这个物理量也会改变。

那你其实就有一个呃partial by party partial t，那这个其实是一个可以说它是它的这个material，derivative的一个欧拉的分量，就是空间中不动，然后随着时间改变而改变。

第二个分量呢是呃，由于你的材料移动了，所以你这个地方你的变化率是什么，变化率是你的u乘上呃，你的这个scale field在空间中，本来关于空间的变化率就是它的gradient对吧。

u点乘上gradient又得到物理量，关于材料呃，由于材料移动所导致的变化，这个material rita有很多名字，然后有的同学喜欢叫他，advertive或者lagrandia或者particle。

他有很多很多的名字，然后我们一般呢这边这个课里面，我们就叫它material，我们不要引入特别多的名字，要不然他会比较confusing，那么常见的例子呢，就比如说温度的变化，温度变化有两种成分呃。

如果是你是把温度计在粒子上面，那你这边啊温度变化，首先它你这个例子不动，它也会变对吧，第一个分量就是怕手t pasht，怕是大t怕是小t，就是他对呃t大t是temperature，小t是时间。

第二个分量呢是以粒子移动导致的变化，那么就是u点乘上gravity对吧，有的时候呢这个material drs也会对速度本身去求，那么如果说你是对速度的x分量需求，那其实你直接套用公式就可以了。

只不过可能有点绕，因为你是呃速度的变化，它本身是老速度的影响，ok刚才讲了这个material derivative，那么我们下面就讲一讲，这个不可压缩的n s方程，这个a商城有很多变种。

这边我们介绍一个比较简单的版本，看起来好像很复杂，但是其实大家不用担心这个方程是什么呢，我们直接看这个公式三，它是什么意思呢，就是说你的流体里面的一个小的粒子，或者一个小一小包材料。

它的速度关于时间的导数有三个成分，第一个成分是由于压强的梯度导致的，这么一个贡献对吧，你如果想如果说你的这个呃一小包材料，它的压强在空间中的变化，那它就会受到压强的这个力，对它的加速导致的影响。

第二个成分是什么，第二个成分它的粘性，这个联系下我们一般在图形学里面是忽略它，因为我们一般希望流体这个效果越不练越好，然后因为越不连以后你就更更为突出各种窝呀，各种看起来非常炫酷的效果。

所以我们一般会把这个联系项给丢了，所以这一项大家其实不用考虑，然后还有一个g g是什么呢，g就是重力加速度，对吧嗯我们再看一下这个方程五，由于我们这个是一个不可压缩的流体，那么不可压缩是什么意思。

不可压缩，就是说你的速度场不能有散度对吧，如果你速度场有散度的话，那你局部就会由于材料的过度流入或者流出，导致这部分的流体被压缩或者延伸，被压缩或者延伸有什么后果，就是它的密度会变化。

但是我们这边又是一个不可压缩的流体，所以它的速度的反应应该是零，那么实际情况下我们怎么来求解，这么呃一个不可压缩n s方程吗，看起来好像很复杂，但是我们一旦把它分开来，分成好几份，它就不再复杂了。

首先我们会做一个叫operator speaking的操作，这个名字又是听起来很复杂，但是大家可以去看一下这个more details，它面会介绍这个operator spleeding是什么意思啊。

它的严格的定义，但是我们这大家只要知道operator stating呢，就是把一个呃，把一个关于时间的pd给它拆成若干个部分，使得每一个部分都非常的简单，这就是它直观的含义。

ok那么我们怎么来split呢，首先我们去把它这个其中的advection这一部分，给它提取出来，也就是说我们把原来这个方程，我把右边右端项给忽略掉，我们认为右转项是零，就得到公式六对吧。

然后我们再把这个g的这一部分给他考，考虑进去，就得到了公示期，然后呢，最后我们再把关于压强的贡献，这一部分给它考虑进去，然后我们还要把呃速度场散度为零，这么个条件也考虑进去，就得到了我们的公式吧。

ok那么我们接下来会分别讲公式六和公式八，公示期非常非常简单，公示期其实往往就是给你的速度上，直接加速一个呃，直接提添加一个重力加速度，然后ok我看到有很多同学说非常非常卡，我想想我有没有什么办法。

我这好像也没有什么办法，我看看我能不能呃，是不是我的wifi的问题呢，应该应该还好，ok我们就继续讲吧，我觉得好像如果卡的同学能不能试一试，换换线路，然后我这就有的同学说正常，有的同学说卡。

那我猜应该我这应该还好，然后ok我应该暂时没有什么能做的，所以就不好意思了，我们就先继续，好，那么我们刚才讲了operator beating，我们刚才讲公式678对吧，我们把原来一个很复杂的一个pd。

给它拆成了好几个部分，那现在呢我们得到了这么1233个部分，外力和projection，那么我们可以再对它进行时间的离散化，那我们就得到了对于每一个time step，我们有三步，第一步呢。

所谓的advection听起来是个很高大上的词，其实含义非常简单，它就是移动我们的流体的长这么呃，或者说呢我们通过当前的速度场求解出来，下一个时间步的速度场的一个初始版本。

但这个初始版本是没有考虑到外力，没有考虑到压强的，我们只考虑速度场的移动，所以这一步叫做action o，那么第二步呢，第二步是在我们刚才求出来的速度上，的基础之上，我们再去加上外力。

那就得到了一个u star star，是两个star，那么第三步我们这个u star star由于经过了，即使他即使u本身是没有散度的，由于经过了the vaction和外力的修改。

他自己可能会产生一些散度对吧，那么我们怎么办呢，我们再进行一步projection，给它加上一个压强的一个作用，这个压强的作用呃，是这个压强当然不是随便抽出来，这个压强得是满足，加上这个压强以后。

你的新的速度场的散度是零，ok那么我们来呃，有了这三步以后，我们得考虑一下，到底怎么样来表示我们的速度场，这边就牵扯到呃我们的一个数据结构，这个叫做实际上非常简单的，就是一个均匀网格。

当然你如果是oa farm模拟流体的话，不一定要网格，你可以用三角网格啊，你可以用这个呃，三维里面你可以用四面体网格，然后你也可以用各种其他的像呃，你也可以用，你如果真的愿意的话，你也可以用例子。

就是例子永远不动，然后你可以用它来build一个欧拉法的方程，那也是可以，但是这些方法都比较复杂，我们今天呢就讲一个非常简单的，这个叫regular gra，均匀网格。

均匀网格呢一般在二维里面就是正方形，三维里面就是六面体，那么如果说我们把空间分成很多个均匀的，正方形的cell，就是cell就是细胞了，我们这边其实是3x3=9个sl对吧，那么这个问题是在于。

我们对于每个cell它自己还是有一个面积的，它里面任意选一点，我们都可以用来存速度的横向分量，纵向分量以及压强，那我们存在哪儿呢，这边就有选择了，最简单的乘法就是全部存在cell的中心。

也就是说偏离横向偏移零点呃，和cl的corner偏离0。5个cel，然后纵向呢也和cell的这个左下角，偏移0。5个cl，那么你把所有东西全部存在中心，如果这样的话呢。

那你去分配这个速度场和压强场的时候，你就可以直接呃在太极里面，其实你直接开这个呃uv p它就是其实ti divr，然后t i。f32 对吧，然后它的shape，其实就和你的cell的大小是一模一样的。

因为你有一个cell，你就有一个压强或者速度的采样，ok然后，稍微更复杂一些的方法呢是用stagger grade，所谓staggrade，就是说我们可以把速度的横向分量，速度的纵向分量以及压强。

存在这个cell的不同的地方，然后这边呢我们红色的存的是速度的横向分量，绿色的点呢存的是速度的纵向分量，蓝色的点存的是压强了，注意这边是写错了，这边不应该是orange，这边应该是blue。

然后如果说你要这么去分配你的速度，场的grade的话，那么你一定要稍微注意一下，因为呃如果说你的网格是3x3，那么你的横向的速度长，它其实是一个4x3的速度，看你的纵向速度长了。

它其实是一个3x4的数组，你看看，因为呃虽然你有三个三个grade，三个每个方向，这里有三个grade，但是你其实如果说存在边上的话，那你其实有四条边对吧，这个我记得小学的时候。

学过一个叫做什么植树问题，还专门研究这种这种这个加1-1的事情，这个可能是童年回忆了，那么压强就不用担心压强还是存在他的cell center的，所以不用担心，那么staggree有什么好处呢。

staggree后面我们算到各种有限差分的时候，它就会有很多的好处，特别是它能够避免一个叫做checkerboard，pattern的这个一个很坏的事情。

如果关于这个checkerboard pattern，大家有兴趣的话，可以去看之前提到的一本书，ok那么我们刚才呃，虽然在每一个cell里面都加了一个采样点，使得我们能知道那一点的速度，那么还有很多。

但毕竟这个excel里面除了一些点以外，还有几乎所有的地方它都是没有数值的，那么怎么办呢，我们怎么样在这个cell里面，所有的地方都把我们，比如说速度横向分量，或者速度纵向分量来产生定义呢。

这个地方就需要用到一个插值技术呃，常见的情况，常见的插值方法是这个双线性插值，by linear interpolation，有的时候有些书上也会讲。

流体里面可能会用by cubic interpolation，但是最近好像不是特别流行，最近好像流行的搞，他就是直接搞一个白linear interpolation。

by linear interpolation，非常容易实现，然后呃在大部分情况下work也挺好的，所以大家一般就用白零点了，什么是blinear interpolation呢。

假设我们在这个一个cel的四个顶点上面，呃，如果我们知道一个一个这个正方形，它四个顶点上的值，比如说呃这个地方还有黄色，绿色，粉红色和蓝色对吧，那么我们如果说需要用这四个顶点上的值。

我们只知道这四个顶点上的值，怎么去重建出这个正方形里面任意一点的值呢，这个其实嗯还是比较直观的，我们其实只需要呃，首先我们只需要把这个点关于嗯，每一个顶点的距离给他算出来对吧。

这个距离是我们是用了一个什么metric，这个距离我们用的metric其实是呃，它关于对角线这个点形成的面积，这个可能稍微有点绕，但是大家可以看底下的这个直观的解释。

我们这个黑点它的数值是怎么加权计算到的呢，它是黄色的点的数值乘上黄色矩形的面积，加上蓝色的点的数值，乘上蓝色矩形的面积，加上绿色点的数字，加上乘上绿色矩形的面积啊，但是这个就是一个简单的加权平均。

其实还是很容易理解，那么做了这一步以后，虽然我们只是在空间里面的很少个离散的点，上面有数值采样，但是我们可以通过双这个双线性插值，得到空间里面几乎所有点的呃，一个数值的近似。

假设我们的物理的这个这个feel，它是一个几乎是连续的，没有这个跳变的话，那这个binary interpolation，还是能够做非常好的工作的，ok那么接下来呢我们就讲一讲advection。

所谓的action啊，就是移动流体的材料，说到reaction呢，这个有非常非常多不同的schemes，然后不同的scheme呢，它其实是在数值粘性，稳定性，性能和复杂程度之间在做不同的取舍。

当然这里面最简单的就是，我们后面第一个讲的叫做similar roger fection，它非常的稳定，但它的数值连性也比较高，什么数值联系呢，我们后面会演示，那是一个稍微复杂一点的方法呢。

叫做micromic或者b f e c c，它是啊3米2个中心的一个升级版本，然后呢最近还有一些新的skin，比如说呃，我们曲子怡同学和张艺兴师兄，提出的这个by mark square。

那也是一个很有意思的陪伴，大家可以去看一看，最近stay up应该是2019的一篇很好的paper，然后还有particle action，就是说我们虽然我们是用欧拉网格，表示背景的速度场。

但是我们可以把很多物理量给它记在例子上面，呃，这个我们在后面讲到，混合欧拉拉格尔说话的时候，我们会提到这个particle action的方法，它其实也有很多很多的变种。

那么今天呢我们就讲一讲最简单的similar，grand vaction和bf c c或者microy，my comic和bf e c c本质上是非常类似的，基本上只有一些系统上的差别。

所以呃大家一般也不做严格的区分，ok那么什么是三bilibrarian roger advx呢，假设呀我们要求下一个时间步x点的速度，那么我们如果知道这个时间不它的呃，上一个时间步的速度场。

我们怎么来求这个下一个时间不x点的速度呢，假设我们的速度场在空间中每一点都是恒定的，并且它从来不变，那么怎么办呢，其实非常非常简单，我们只需要把x点的位置沿着时间，沿着当前的速度上去嗯。

向后推移一个吊塔t，向后推一个吊塔t，这样我们就得到了x这一个点，在上一个时间步的位置对吧，那么由于x啊，它是你可以认为它是一个拉格朗日的一个例子，那么由于我知道他上一个时间步在这个地方。

下一个时间步跑到了这个地方，那现在问题来了，这一个时间不x的速度是多少呢，这个其实非常简单，其实你只要在上一个时间不嗯，在x这一点做一次双线性插值，你就知道了，这个例子在上一个时间步的位置啊。

在上一个时间步的速度对吧，由于它这个粒子的速度，你假设它是不变的，所以在下一个时间步的速度呢，就等于呃这个例子在上一个时间步的位置，在速度原来速度场里面的采样，ok那这个是要实现起来。

其实非常非常的容易，我们对于，如果说我们要对x这个物理量去求similar，graction，把结果存在lex里面，那么怎么办，我们就去loop over x里面所有的下标，然后把new x。

他这个下标的值设成blinear interpolation啊，我设成x值就是x是旧的物理量的值，在这边，在它的历史上的位置去做一个双线性插值，这边呢我们要做一个back trace函数。

back trace函数是什么意思，就是说我们把时间倒着播放，倒着播放，我们就能回到历史，就可以求出x这一个点在上一个实验部的位置，然后我们在上一个实验部去求x的值，就可以得到呃，这个时间播的物理量。

ok可能有点绕，但是呃应该直观，应该直观上还是比较容易理解的，但是呃这个就会导致很多问题对吧，我们由于我们的速度场它并不是一个恒定的值，它会经常变，那你实际上x它在速度场里面的，这个x这个支点。

还有速度场里面移动，并不一定是沿着一条直线，它可能是沿着一条曲线，那么实际上你的这个material parcel，它的移动路径可能非常非常的复杂，你可能会得到一个偏离后的偏离后的结果。

比如说呃你如果这次这个支点x是在这，那么你会计算你获取会计算出它上一个时间步，它的位置是在红色的点，这但是呢实际上他可能是在绿色点缀，那你这个就设置好礼仪，没有1000里了，那么怎么办呢。

队友在讲这个之前，我们先来演示一下。

这个这个误差会产生一个什么样的效果哦，对其实代码现在已经在github上面已经有了，应该是在这，大家可以去论坛里面去点击这个链接，或者你直接到games 201的这个release里面。

可以下载到我刚才演示的这几个几个代码。

ok那么我们来演示一下，你如果看到我的vx里面其实有一些选项，我先把这边来稍微调一下一。

ok那么我们一开始有一个太极的图标，然后我们加一个旋转的速度，成如果说啊我们用一个比较简单的sevlog，range and action的话，如果说我们只用一个非常简单的trace back。

的一个scheme，那就会导致这个素质场，他的你看他这个太极图标就越来越小，越来越小了，这个就是我们的误差导致的。

由于说呃你每次trace spect的时候，你的这个路径都不准确，你其实会得到一个更远一点的值，导致你的这个每次做完一留以后，做完advaction以后，你的图标就越来越小，当然还有一个问题是什么。

它这个图标会越来越模糊对吧，这个越来越模糊的问题，我们后面再解决，我们先来解决一下，这个图标变得越来越小的问题，如果说你呃如果说你有一坨流体，它就是这么顺时针的旋转的话，那么你如果说我们这个太极图标呢。

就是表示一个速度场，白色表示高温度呃，黑色表示低的温度，我们假设速度场不会defuse，那么我们呃呃我刚才说错了，应该是温度场，白色是高温度，然后黑色是低温度，我们假设这个温度场是不会defuse的话。

呃，那么我们当然不希望这个温度场，随着你流体的旋转变得越来越小，这个就会就实际上是数值误差对吧，那怎么办呢，我们其实有更好的backtrace的方法。

这个其实就变成一个initial value problem，它是一个数值o d e，学过数值分析的同学，应该都知道这个该怎么处理，我们刚才啊其实演示的时候，我得先把这个背景的程序给关了，要不然太卡了。

ok我们刚才演示的时候呢，其实用的是forward order，然后它其实很简单，我们要求这个p在上一个时间步的位置，怎么办呢，我们把p减掉，p在当前时间步的位置乘上dt对吧，那么我们其实就回到过去了。

那么这个叫做fooiler，他其实呃注意这边，其实这个负号是因为我们在往后面往历史去推，而不是往未来推，如果是往未来推的话，这个应该是一个，正好这边里面所有的符号都是，由于我们在going back。

我们在重新播放，我们的返回播放我们的物理模拟导致的，那么怎么得得到一个更精确的一个estimation呢，其实直观上也非常容易，我们可以用ark to game或者叫显示终点法，显示充电法干啥呢。

我们先把这个p呃往回回播放，0。5个dt，我们就得到一个两帧之间的中间点，然后呢，我们在这个中间点这个地方采样一个速度，然后用这个速度去play back啊，t的位置，ok那么其实你可以想象。

如果用了终点法的话，我们的对于速度的第二次，对于速度的estimation就会更准确一点，因为呃你如果一个区间你要近似它，你怎么近似，可以用开头的，可以用结尾，但是最好的还是用中间的。

如果你这个区间跟上面的数值，是一个线性函数的话对吧。

那么我们先把它换成我们试一试，大家可以改一改这个scheme呃，我刚才rk是一，我们现在把它换成啊换成二，这个注释应该在这换成二呢，我们机器实在用ark to cheme。

你可以看到我这边有个叫backtrace的函数，然后它里面有一坨ta。static对吧，然后ti。sestatic，其实就是在编译期做一些判断，那么如果当rp是q的话。

我们就是在用中点法进行呃向后的追踪，向历史的追踪。

我看看重点法效果怎么样，ok你可以看到用户终解法以后，我们基本上没有了，刚才那种越转越小的问题对吧，虽然它还是越来越糊，但是越来越糊，这个是另外一个问题，我们先把这个越转越小的问题解决。

这是2k two，一般来说2k two呢有的时候已经够了。

我们也可以试试2k three。

这是2k three，基本上和2k two没有什么区别。

ok然后，2k three是干啥的，2k4 其实就是一个升级版的2k two，它更复杂一点，它首先呢呃往回走0。5dt得到一个p1 ，然后再p处呢这边再sample一次，然后得到一个呃，得到一个v2 。

然后用这个v2 呢再去play back一次，得到一个p2 ，然后再再三步一次得到一个v3 ，然后你最后真实的速度的estimation，其实就是啊，把之前你得到这些速度进行一个加权平均。

然后这个全职是不是随便写的，这个全职是有技巧的呃，它有一个推导的方法，这个大家可以自己去研究，但实际情况中一般2k q呃，一般就差不多了，有的时候用2k4 台，很少情况下必须要用2k4 。

但是随着你的这个格式越来越复杂，你的计算量也会变得越来越大，对吧呃这边就也是有一个取舍了，然后我们刚才讲到啊，如果用sama rgx，它会越来越糊，为什么会越来越糊呢。

还记得我们之前讲过blinder interpolation吗，由于我们在做interpolation的时候，我们会做甲醛平均，那平均了以后，其实本来你如果有一个很sharp的一个变化，你做了平均以后。

它就变糊了对吧，你本来如果这个端点左边是零，右边是一，你取在中间取一个，你就得到0。5，那你这个速度场就会很自然的变糊，那么怎么样解决这个变糊的问题呢，这边就要用到一个稍微高级一点的呃格式。

这个叫做bf c c或者叫做micromic，然后他们两个刚才说的其实非常的像，然后大家很很多文献里面其实会混淆使用，他们俩呃，这边我们就讲的其实是底下代码其实是microy，然后啊。

他你可认为他的hilid其实就是b f c c，什么是b f e c c呢，b b c它的全名叫做back and forth arrow，compensation and correction。

什么意思呢，其实想法非常简单，你可以看到它其实调用了这个sl是三倍，lol running action，其实x新就是我们刚才的sema grand，ok，那么我们用x星啊，先回到过去求到一个呃。

把x做一次of action得到是得到一个呃x x，然后呢我们再反过来做一次，我们再把x音做一次ction，当然这边我传给他的是负的dt，就得到一个x星，x星星对吧，那么理论上来说。

如果说我们这个semi oraction，它是没有任何误差的话，那么我们把这个速度长先往呃往前推一个dt，再往后推一个dt，我们是不是应该得到x星和x星，星和x它应该是一样的对吧。

因为呃如果说你的action是完美的话，那你其实往前推dt，再往后推推dt，你应该是能保持原样的，但是实际上它不一样，那不一样以后呢，为什么会不一样呢，是因为advection error。

那么我们把它俩做差，然后再取个1/2，就可以得到我们对advection error一个估计，ok我们我们得到这个reaction error估计以后呢，得到一个x final x final。

其实就是我们的最后的速度长，ok但这边你要小心一点，就是说如果说呃你这个correction本身啊，可能会导致你的the vaction scheme产生overshooting。

这个会产生一些artifacts，然后我们现在来看一看bf e c c的效果。

我们大家如果下载了下载了这个代码，可以进行一下修改，比如说我们可以把这个u s m c换成q对吧。

mc就是micromic，然后可以看到我们一旦用了micromascheme以后呀，这个基本上就没有这种变糊的问题了，基本上没有变糊的问题，你可以看到我，我现在按一下空格可以暂停，你可以看到呢。

它这个圆这边部分这个部分基本上没有变弧，它运动的部分还是稍微有点变化，但这个和三倍lab软件比起来已经好了很多很多。

我们来试一试，再看一看3米的高中卷，可以做一下对比。

我刚才把uc mc设成了false，然后现在是在用sam两个range，你可以看到，特别是这个圆的这个边界啊。

它已经变得非常的糊了。

在上转一会儿，他直接就，你让他转一会儿，这个地方其实已经不太看得清了，这个不太看得清以后会有什么后果呢，一旦不但看呃，这个是一个其实是一个颜色，一个黑白的一个颜色场。

如果说大家是在速度场上面做这个操作的话，它一旦变模糊，就说明你产生了数值的联系对吧，数值的联系就是速度场呃，速度场在模拟的过程中，它会产生一些nero的一些diffusion的这个问题。

它有产生这个这个diffusion以后，它就会导致能量的衰减，量衰减以后，这个流体啊就没有之前看起来那么energetic，所以就会有问题，就会导致你的流体看起来非常黏，你看转一会以后。

他这个就基本上不太看得清。

刚才讲到accomic里面有一个overshooting的问题。

我们来看一看，我们还是现在有my home，这个要有一个clipping，keeping是干啥呢，keeping是我们做完这个呃micro mac以后，我们需要去防止它产生一些非常大的值。

为什么会产生非常大的值呢，因为它这个地方它会根据你的呃这个reaction error，让你在边界上产生一些呃，像吉布斯现象那样的问题，然后如果你不去做这个一个呃guarding的话。

它就会导致产生一些artifact，我们来看一看吧，如果说我们把这个clipping给他去了，你可以看到呢，它产生了一些很看起来不是令人满意的，这个波纹状的这些artifact啊。

所以大家记得如果写macrome的话，要记得去做一下clipping。

但这个代码都在ction点拍里面有，大家可以下载下来，自己玩一玩，ok那么刚才就是direction，那么接下来呢我们讲一讲projection，projection啊，其实它非常简单。

projection是主要要做什么事呢，你得去求出一个合适的标量场p，使得这个标量场p的梯度呃，施加在这个u上面以后，你的u变成了一个没有散度的一个速度，成什么意思呢。

我们用有限差分在空间上面展开一下啊，我们假设投影后的这个速度场，project以后的速度场叫做u star，那么u star减去u这边应该是除以个delta t，等于什么，等于负的肉分之一啊。

gravt对吧，我把这个交叉t呢从u2 减去u下面，移到了这边来，呃当然我们要满足usr的导数啊，usd这个散注是零对吧，那么我们可以做一做一项得到了公式十，然后我们在对公式时的左半部分。

这一部分两边同时求散度，为什么两边同时求散度呢，因为求了散度以后，我们左边就得到了一个这个呃u s2 的散度，u s散度是零对吧，我们希望把这个usr给他消了，同时写了散了以后，我们知道左边变成零。

右边呢呃散出operator是一个linuerdo，可以把它塞进去，得到了一个u的散度，减去掉他t除以肉，然后然后有一个这么一个算符，这个是一个laplacial period，后面我们会提到它。

其实就是说对于p的gradient求一个散度，ok那我们再做一个最后的一项，得到了公式13，公式13呢其实就是呃，p的梯度的散度等于肉除以调tt，然后再乘上u的散度，那么我们需要求出一个p。

它满足公式13，这个公式公式13啊，到这变成公式14，然后公式14，它其实呃是一个非常非常常见的一个形式，它叫做pssecreation，它呃你如果把右端项做成f的话。

那你经常可以把它写成这个delta p等于f，但是这个da大家注意，这个delta dela p和delta t dotx不一样，它不是说p的差值，它是个da，它是laptoperator。

然后他有的时候呢会写作lab平方，有的时候会把它写成这个呃，中间一个点乘的一个形式，两个倒三角，这门课里啊，由于我们用deltt deltx来表示呃，差值表示这个difference。

所以我们就不用这个调查符号来表示laptoperator，这样避免一些歧义，如果说呃f等于零呢，这个方程又被称为拉普拉斯方程，但是我们这个课里面f一般是不是零分，因为我们一般是用这个泊松方程。

是用来求解速入场的压强来投影这个速度场，所以我们的一般这个f啊它是有良好的，f其实就是公式14的右半部分，ok那么我们下面下面要做一个事情，还记得我们的物理量，都是存在一个背景的网格上面的吗。

那么我们要建立一个离散的方程，来去近似这么个离呃连续的方程，ok那么我们要建立一个线性方程组，这个线性方程组呢，它满足a乘上p等于b的形式，那么a和b分别是怎么定义的呢，ok那既然我们的方程式公式16。

那么我们看一看，我们首先我们看一看怎么来近似这个p的拉plus，operator呃，这个拉普拉斯p在ig点的求职，这个我们后面后面一个slides会讲到，为什么会是这么一坨复杂的公式。

大家这边只要记住它是用p周围的四个点的值，包括他自己就是一共是五个点，来求出p在这一点的拉普拉斯freer的一个近似，这个叫做five points dl，号为的five points dl。

然后还有一个要求的是什么呢，还有一个要求的是u的散度对吧，u的散度怎么求呢，u的散度其实也可以用啊，u周围的四个采样点的值来算出来，这个我们下一页slides会介绍，为什么是这么个公式。

大家看这个信息系统啊，这个线性度还是挺大的对吧，如果说你的呃网格是n乘以m的话，那你这个线性系统，由于你的p p它是有n乘m个采样点，你的b也有n乘m采样点，那你的a呢其实就是n乘m再乘上n乘m。

如果你是三维的话呢，你的a这个矩阵就是n乘m乘k，再乘上n乘m乘k这么个元素，它其实就是一个非常非常大的矩阵，好在它是稀疏的，我们后面会讲到怎么样去呃，求这个矩阵a，ok那么我们先来看一看。

为什么u的散度会这么被近似呢，哦注意这边少了一个e除以dtx，这个u的散度是什么意思呢，直观上理解就是说这个cell里面呃，速度的散度意味着什么，意味着你这个流入或者流出的流体的量对吧。

如果你是流出的话，那你的这个散度为正，如果你是流入的话，你散度为负，ok那么我们怎么样来求出呃这一个点，它流出流体的速率呢，还记得我们之前讲了staggrid吗，如果说我们用了staggrid。

那么我们的速度的横向的分量，是记在这个cell的两条竖边上面，速度的纵向分量呢，是记在这个cell的两条横边上面，ok那么如果说我们知道这个这一条边上，它的速度是呃是uxi加1j的。

这个它的流体是以这个速度往外面流，那么它这个cell它这个这条边对于整个cell里面，流体流入或者流出的贡献，其实就是吊塔x乘上ux a加e j，那么我们还得把它减去它左边的这个部分对吧。

你左边是流进来的，这边是流出去的，流进来的得减掉，流出去的得加上，那么纵向的两条边啊，这个横向两条边其实也是一样，它的纵方向的流体对于流体速度，对于这个cell的它的嗯密度的贡献。

他其实是首先你得加上他这个流出去的对吧，然后你再把它减掉，它流进来的，其实这么1+1减，不好意思你你就得到了他的呃，速度场的散度，那么压强大plus of operator是怎么算的呢，其实也非常类似。

如果说我们要求在style中间，这一个点的压强的lapoperator，那么用类似的推导方法，我们也可以先把这一个cell，它的四条边上面压强的梯度给他算出来对吧，这一个点它的压强梯度是什么呢。

应该是上面这个点减去下面这个点，除上吊塔x就得到这个点的梯度呃，这条边上面的梯度类似呢，你四条边的梯度都可以算出来，然后你要求它散路对吧，你要求拉散路怎么办呢，就和上面一张size一模一样。

那么你做了这个操作，再进行一点点化简，你就可以得到这个嗯p的laptoperator，然后-4 p i g，然后加上它周围的这些cell的值，那么现在问题来了。

如果说我们这个流体假设它右边这三个cell是啊，是一面是一个固体，假设我们的流体就左边两列，第三列呢是一个固体，那么这边我们就得对它进行修改，这个叫做呃这个就是它的边界条件的问题。

如果说他右边这三个cell呢是空气的话，我们就叫做它是一个开放边界条件，开放边界条件其实非常容易，我们就可以把p认为在这些cl里面是零，那么在下面这个方程里面，我们只要把t的对应的值换成零就可以了。

那么如果是它是一个固体固体的话，也就是说固体什么意思呢，也就是说你这一面它的压强差或者这一面呢，它的速度一定是零对吧，那他的速度永远是零了，那你的这个压强差呃，就其实压强的梯度应该也是永远失灵。

那么如果说他有一面，你会发现它呃，如果你发现它有一面和一个固体相邻，那你其实首先你得把这一面的这个，p i g的贡献给去了，那除了去了他呢，你还得把这个四给它减掉，一变成三对吧。

因为你中间的这个它的贡献，也因为你的右边这条边是一个固体，而去掉这个推导，大家可以去看之前提到的一本书，那本书里面robert ron，那本书里面讲的其实非常的细，这边由于时间关系，我就不继续深入了。

ok刚才讲了怎么去建立一个方程，a乘上p等于b对吧，那么这个方程我刚才提到，他是非常非常large scale的一个线性方程组，那么我们这怎么去检查呢，首先啊，我们来看一看这个20世纪最重要的十个算法。

但这个是仁者见仁，智者见智，这个只是呃有一个期刊评，评选出来这么十个算法，十个算法我觉得还是比较有意思的，你可以看到呀，我这边加粗了其中的四个算法，这四个算法有什么特点呢，他们都可以用来解线性方程组呃。

我们今天要讲这个quite of subspace，iteration method，然后它其实是历史非常非常悠久的一个算法了，然后呢除了用iteration method解linear system。

你还可以去做matrix decomposition对吧，你可以去prefect orize你的你的矩阵，然后呢你还可以用，你如果是要解泊松方程了，其实你还可以用快速布列变换f t。

然后呢你还可以用fast multiple去做preconditioner，所以其实你看这些算法里面，其实很多都是和线性方程组有关的，就是ax等于b这么个线性方程组，它背后有这么多的学问，嗯顺便说一句。

这个simplex methods for india programing，其实也可以用来解rage body simulation，然后my chocolate originary。

mtchery m c m c其实也可以用来做衰做渲染，其实图形学里面用到这些算法还是很多的，另外一方面也确实说明，这个算法确实在各个领域都有很好的应用，ok嗯对，其实由于我自己把虎牙的直播给关了。

所以我现在其实看不到大家弹幕里面的问题，如果大家有问题呢，请在课程的微信群里面提出，这样我就可以或者在games的群里面提出，这样我就可以看到了，我不知道大家现在暂时有没有什么问题，因为我把网页关了。

我自己也看不到，所以大家可以发到微信群里面，这样我就知道了，ok那么回到我们解信息系统，很多物理引擎啊，最后到最后都都转化为一个解，一个很大信息系统的问题，ax等于b然后怎么去检查呢，有很多种解法。

如果你线性系统不是特别大，你可以用partil这个direct solver呃，其中呢最有名的就是partil，然后大家也非常喜欢用它，毕竟你只要把这个信息系统给他的，把这个矩阵的气候形式给他构造出来。

然后喂给parties，把这就给你很稳定的给它解出来，然后所以大家很多时候会很偏爱，这么这么个方法，因为它省事，而且性能也不低，其实在系统比较小的时候，不是特别大的时候还是很管用。

除了direct server呢，还有一种一类server叫做iterative ser，direct ver呢会一下子给你把解算出来，iterative solver呢它会呃得到一个解。

然后在这个解上面，然后呢这iteriterative ser有好几种呃，很多种变种，比如说你可以有gala，还有jacob，然后这两个是非常常用的，除了他们还有什么s o r i。

s o r之类的各种呃，其他的itering server，其中呢damage隔壁我们之前已经讲过了，今天呢我们主要讲的是这个precondition country。

country regredients，然后它是快乐subspace over其中的一种，这边要提到呀，一个好的emicrosover，通常呢是有很多不同的solver给组合成的。

比如说图形学里面要求这个泊松方程，一般会怎么做呢，大家外面会用conjugate gradients，然后用multies or preconditioner。

让multigrade的每一层呢用djacoby domoving，然后呢在multivate最下面一层，用partial这个direct sober去进行direct off。

这个基本上已经是图形学里面，求泊松方程的标配了，呃这个里面有很多名词，我后面会提到，这边大家只要记住它是呃，solver之间可以进行组合，那么现在问题来了，我们怎么去表示这个a呢，刚才我提到呀。

如果在三维free simulation里面，这个a可能会非常非常大对吧，因为你每一个vocal都得在a里面占一行，那么但是a有一个通常来说，如果你运气好的话，它有一些好的性质。

比如说呢它通常是sparse，如果运气好呢，它还是对称的，并且是镇定，叫做对称正定或者叫s p d，那么你怎么去存a呢，a说白了其实就是个二维数组对吧，你最简单的方式你就把它存成一个稠密的矩阵。

一个dance matrix，这个通常在a比较小的时候是work的，但是它不是特别scalable，比如说你如果这边是1024，是1024，你还可以存，如果说你有一兆个呃voxel。

一兆或者一兆个pixel，那你就存不下了，你都没有空间存它，这个时候呢往往就要用到一个sparse，matrix的一个表示，用sparse matrix来存的a，为什么呢，因为虽然你a的每一行。

可能比如说你有一兆个元素，但其中呢很可能只有十个元素是非零的，其他全部都是零对吧，你如果你全部都是零，那你就不用存它了，然后sparse matrix有很多很多种存档格式啊，什么csr啊，coo啊。

然后这边就不具体展开了，然后在图形学里面一般的常见的搞法呢，我们是不存这个a很不存这个a什么意思，这个我们后面会提到这个叫做一个matrix，three的搞法，通常呢这个是最高效的方法，对于存储来说。

由于不需要存a的entry了，所以它能省很多内存，对于计算来说，大家不要觉得计算是加减乘除公式计算，你从memory fetch data这个也是计算一部分，而且在现代计算机体系结构里面啊。

memory bandwidth和flops比起来是非常昂贵的，一个一个这个操作，而是非常昂贵昂贵的资源，因为你内存很慢，cpu很快，这个我们后面第九讲会提到高性能物理模拟。

这个时候你就要考虑到我到底是在线的，把这个值算出来，还是说嗯我去把这个值预先算好，放在内存里面，然后每次去取这两个，是这两个设计在不同的情况下会各有优劣，但是对于呃流体模拟里面的泊松方程。

一般来说大家都是不存的，直接用magic three的方法去表示这个a，ok那么具体matrix free是什么意思呢，我们后面会提到，这边我们讲一讲clot subsidivers。

然后quite a subsisubsive solvers，基本上是最有效的linear systems office之一了，然后它其中又有很多很多个变形。

那最知名的变形叫做conjugate gradient，共轭梯度法，或者一般会被缩写缩写成cg跟computer graphics，其实呃那个cg不是一个cg，然后在graphics里面呢。

除了cg以外，大家还会用，有的时候会用cr，还会有这个gm s，然后可能很少很少情况会用搬到bcd step，大家这边呢一开始的话只要掌握cg就可以，因为cg确实在很多情况下已经足够了。

如果大家想学更多的关于quite of super，super space over的话，有一个书呢叫做introduction to conjugate。

gredimethod with all day啊，就这个不过梯度无痛教程，这个还是一个写的非常好的书，这个书已经出来呃，几十年了，然后呃这个应该是免费，这个应该是一个免费的书。

因为他直接挂在自己的主页上了，大家可以下载下来看一看啊。

虽然他自己声称是一个无痛教程，你如果真的要把它啃下来，还是要费点时间，我自己之前看的时候，说实话也费了一些时间，并不并不像他说的那样完全是无痛，但是他标题也没有标题党，他说这个是没有剧透的学习。

你这个只是呃学习嘛，总是要有一点痛苦的，他只是把你痛苦减轻一点，这是一个非常写的非常好的书，大家可以去看一看，关于什么是快乐subspace，我这边就不继续去深入了。

然后对于conjugate gradients呢，其实它实现起来非常的容易，其实也就这么个十几行代码，ok然后如果说你一开始没有办法理解，conjugate gradients在做什么，完全没有关系。

你把这个十几行代码啊记好或者给他抄好，那你基本上你soba就是work，然后呢你当你得到结果以后，你再去慢慢研究它，为什么是这么个结果就可以了啊，这个确实一开始要掌握的话，还是需要一些理解的。

ok那么这边大家可以看到呃，他我大概讲一讲这个算法在干啥吧，它其实就是一开始先算一个residual，然后这个p呢这个p和刚才压强不一样，这个p呃是我们的沿着这个解的这个search的方向。

然后x呢是我们对当前的解的一个estimation，然后我们每次呢这个有一个循环对吧，我们每一次在这个循环里面，我们去estimate一下x应该加上多少分的p，使得x能够更接近我们的解。

然后呃我们更新了x以后，这个r也得更新对吧，你看这个r的定义r等于b减去a x，所以更新了x呢，我们r也得更新，而更新了以后呢，如果r已经足够小了，我们就中断。

然后呢我们在estimate一个叫做beta，然后用这个贝塔p和和r去更新新的p，得到下一步的这个p的这个search direction，然后到下一个循环里面，你看我们用这个p来更新x了。

所以这个是一个不断的迭代的一个过程，当然中间这些为什么阿尔法这么取，为什么贝塔这么取，这个就需要一些比较深入的理解，直观上来讲，他就是在做一个能量最小化。

然后他还是选择一个在一个快乐suspace里面去做，能量最小化，还是一个很有意思的一个过程，大家可以去看那本书，那本书里面讲的非常的清楚，但是那本书作者自己也很坦然，他说他写了本书的动机。

就是他自己一开始看了很多共轭梯度教程，一个都不懂，然后他最后最后就自己写出了一个教程啊，所以这跟我记录大家一开始不能理解，我不用担心，大家都是都是不能理解的。

这个本身就是一个很computing的一个过程，o，说到迭代法求解线性方程组，iterative method，他的convergence rate和一个指标非常联系密切，这个指标是什么呢。

是condition number，什么是condition number呢，这个就涉及到矩阵的ion system，我们来回顾一下什么是特征值与特征向量呀，如果说我们有这么个非零向量x。

使得ax等于lambda x，也就是说呢a作用在这个x上面，会把这个x方向不改变，而长度缩呃，缩放拉姆达倍，那么lambda就是a的一个特征值，x就是a的一个特征向量，一个矩阵啊。

可能有很多很多个特征值对吧，然后呢呃condition number是什么意思呢，condition number就是这个矩阵的最大的特征值，除掉最小特征值，这个直观上的理解。

就是说这个矩阵如果作用到一个向量上面，它如果有一些方向把这个向量拉的非常长，有一些方向把又把这个向量缩的非常短，你就会觉得这个矩阵看起来非常的厉害，他的condition number就会非常大。

当然这个是直观的理解，然后要注意的是，condition number有很多很多种定义，我们这边是采取了，通常对s p d矩阵的一个定义，你可能在其他地方呢会看到其他的定义。

但总体来说呢肯定是number越小，它收敛得越快，另外一个呃，能够提高你的这个刚才我们讲的肯定是number，越小才说的越快，那么另外一个能够提高italy solver的一个。

收敛速度的一个trick叫做什么，叫做one starting，那这个是什么概念呢，如果说解这个literally dissolver的时候，你的initial gas已经和你真正的解，非常接近的话。

这个叫做warm starting，在流体模拟里面的warm starting是什么呃，大概是怎么做呢，就是说你这个linear system，解的是每一步的压强场对吧，那么我们有一个sumption。

就是说你上一部的压强场，很可能和这一步的压箱厂差别不是特别大，那咋办，我们把上一步的压箱厂拿过来，当做这一步的压箱厂的初始初始的猜测，然后进行迭代，这样能够显著的在某些情况下能够显著的减少。

你需要迭代收敛的时候需要的迭代次数，当然实际上情况我自己的测试发现，这个trick对于jacob啊，搞赛道啊，或者country gradle是work work，非常非常好的。

但是对于multigrade就word不是特别好，具体的我没有研究过，我觉得直观上的理解可能是multigrade，他的preconditioner distort，你的你的这个空间。

这边我们来演示一下，这边是一个free sober对吧，这是我很早之前写的，你可以看到这边有一个argument叫做warm starting，如果说我们把warm starting给它去了。

你就会看到我这边当然没有迭代到熟练了，我这边是呃，你看我这个有jcopy iteration，我只是用了固定次数的joperation，如果说我把它调到34，你可以看到它的体积会不断的缩小。

为什么体积会缩小，体积会缩小，就是因为你的压强场没有球好，如果你的压强场求好的话，那你经过压强这个投影以后，你的速度场是五散的，速度场无散以后，你的体积就不会放大或者缩小对吧。

但如果你压向上球的不好的话，那你的速度上就会放大或者缩小，那么如果说我把warm starting开下来，我就需要很少的迭代次数啊，我就可以把压箱乘求了呃，迭代到一个非常准确的一个解。

这样的话呢嗯我的这个素质场就会更加的不散，效果就会更好，那么刚才讲到preconditioning啊，刚才讲到这个condition number有没有什么办法呢。

我们可以让a的condition number小一点，其实你让a的肯定是number变小了，那可能是没有办法，因为你不能变a，我是不是又把这个为什么这个风扇这么长，ok，不好意思。

我这个电脑它一跑gpu的程序，它就会风扇就会转起来，大家可能会听到一些噪声，但是我声音大一点吧，刚才讲到呃，a的condition number越大，就会导致迭代法需要更长的时间，或者更多的迭代次数。

有没有什么办法，我们能让a的肯定是那么变小一点的，这个叫做这个技巧，叫做preconditionally，注意啊，我们呃去solve ax等于b的时候对，如果说m是一个可逆矩阵。

那么solve ax等于b，其实等价于我们sob一个m-1，就是m的力乘上a，再乘上x等于m的逆乘上b对吧，我把ax等于b左右各乘上m的力，这个为什么会有帮助呢，这个因为有的时候md乘上a。

会有一个更小的这个condition number呃，和a api，它的可能性当中可能更小，它收敛就更快，或者说呢嗯他的icon value classroom比a自己更更好。

这个icon value classroom，这个是一个country regredient soga，里面决定它收敛迭代次数的一个metric，这个大家可以去看之前的那本cg的书，可这边有个问题了。

有同学会问，为什么我们不直接让m等于a呢，啊这个其实就是呃很有意思的一个问题，就因为你如果那就是m等于a的话，那你其实你这个方程其实就变成了，x等于md乘上b对吧，那其实就不用解了。

但为什么我们不能让m等于a，还是因为一旦em等，如果你都有m等于a了，那你就能invert这个a了，那你就不需要这个preconditioner对吧，这个就是呃就是因为我们没有办法算出m的。

我们就是因为我们没有办法算出a的逆，所以我们才需要做preconditioning，如果我们能让m等于a，那我们就能算出m d m d m等于a，那么m就是a的，你就是因为md算不出来。

就是因为a的d算不出来，我们才需要precondition 0，所以啊一般来说m是取一个近似的，和a近似的一个矩阵，而且他很容易求利o，那么刚才提到m一般是取一个和a比较接近，又比较容易求逆的矩阵。

哪些矩阵比较容易求力呢，对角阵很容易求力对吧，对角阵你直接把他对角线给他求个力，就得到了这个取证里，所以呢这个就引入了这个就造成了一个pregnition，它叫做jacob pregnition。

或者有的地方叫做diagonal preconditioner，还有的地方会用incompleto ley，然后呃这也是一个常见的一个precondition。

这个大家可以去看robert bryson的一本流体的书，里面会讲到sap，cg就是用incomplete lesky来做padditioning，今天我们主要讲一讲这个multique。

因为multique现在基本上已经成为同心学里面，求压力投影的这个标配了，大家还是要能够实现mt龟，猫龟还是很容易实现，而且效果非常非常的好，然后大家不要觉得preconditioner。

一定得能写成一个矩阵的形式，你如果有一个fast multiple，你也可以把它用作pregnitier，他只要是一个general linear preter就可以了，那么我们来讲一讲mi rs。

刚才讲到了mti que，它是一个收敛非常快的一个preconditioner，他其实近似的去invert了a，为什么他能够近似invert a，而且又算得非常快呢。

因为大家想如果说你去解建你建立的这个方程，如果说呃你你只在原来的网格上建立这个方程，那你的这个还记得我们对于laplus，operator的这个discretization嘛。

我们是关于一个cel和它周围的四个cl呃，做的这么一个dispensation，那么他这个每一次做一次这个零点operator，他能看到的cell是非常小，只有五个cell对吧。

但是很多时候你会发现我们的这个呃，速度的散度，它可能是一个非常large scale，非常low frequency的一个速度散步，这个low frequency的这个速度散度。

就没有办法被这个一个senso给它capture了，那么怎么办，我们可以搞一个multi resolution的一个方法，我们把这个recedure这个这个这个余量给它。

不断的做着做一次restriction操作，什么意思呢，我们把四个格子给它合并成一个格子，合并成一个格子呢，然后我们在这个格子上面再做gsmoothing。

然后呢我们做在computers residual，然后再respect，然后最后我们建立了一系列的网格，从最细的网格到最初的网格，到了最初的网格呢以后呢，我们就可以这个一般信息系统就非常小了。

小的信息图怎么算比较快，之前我们还提到对吧，一般就是用direct solver或者用或者用potil，然后当然如果你懒的话，你直接在小的信息系统上面，直接用jcb叠加个240，那也没有问题。

那么这个最对highliywood idea，其实就是说我们把这个residual grade，给他做一次motor solution的处理，然后使得我们可以在多个尺度上面去啊。

消灭我们的residual，这个multi grade它有很多很多的设计空间，它呃刚才看到了一张图，大家也会觉得非常复杂对吧，你这个有有restriction，有prolongation。

有这个v cycle w cf，cyle有好多种不同的递归的方式，这个大家可以呃，我后面会推荐一本书，大家可以去看一看，然后smoother呢刚才我提到好的linear system solver。

他很多时候是由比较简单的linear system server，给它compose起来的，那你smooth里面可以用gf，可用jacoby，可以用dap，jacob和s都可以。

可以是各种各样的各种各样的smoother，然后呢，还有你可能会有很多个很多不同的层对吧，你的层你这个到底选多少层比较好呢，如果你层数太少，那你这个对cos这个level，他的vox数目会非常非常多。

那解解会很慢，如果你存太多了以后，他做了cos以后，那它就会逐渐的失去了原来几个星期的把握，所以层太多了也不好，这个就有一个trade off了，一般来说大家的搞法就是说不断做cosine。

直到你的最cos的一层the box数目，能够被direct solver高效的soft的时候，你就不再继续做cos 0，然后最底层的solver用什么，可以用bf force riobe。

我刚才说的偷懒的方法，或者使用disober，然后呢刚才讲到这个smoothing iteration，到底每一层啊拍了多少个，这又是一个需要设计的地方。

然后呢还有cos和boundary condition，怎么去handle guerin cos啊，sei break mtv的呀，这个后面的话题就很多很多很多了，我们就不展开了，当然说这么多。

其实你要实现起来并不并不难啊，我们有一个参考的实现，一般来说呢我刚才提到，如果你要解personequation，如果你是在regular grade上面来做，那基本上motique是少不了了。

现在你要发一个c graph，你如果要接persn equation，你不用multi gram的，别人都会质疑你这个性能是不是有问题，呃这边呢推荐两个材料。

一个一个就是multique tutorial，这个是非常好的一个教程，然后呢另外一个是graphics里面啊，这个应该不是最早的multigrade paper，但是它非常容易understand。

它就是一个，其实太极里面也是有multi grade的例子的，我们来看一看这个呃，mtv很少单独用来做sober，一般大家都是用来把它用来做cg的preconditioner。

所以最后就会得到一个东西叫做multi grade。

precondition and company regredients mg pcg，我们来跑一跑这个m tp c接，诶advance，这个声音应该是一个三维的播送方程，server。

这个每个iteration啊，就是一个共轭梯度的iteration，你看到这mtv的产生了很多的kernel嗯，我们来看看他刚才的输出啊，在这儿还不太会找到，如果用了mg pcg。

我们只需要六个iteration，就可以把residual从这么大一个数，缩到这么小的一个数，当然m g p cg，你用m g p c g以后，你的每个intuition也会变得很贵很多。

但是总体来说呃，减少的itation的数量还是能够抵消到，你每个h是变贵的成本，所以大家基本上解psy equation都是用mg p c g，呃，我其实可以modify一下mg pcg。

我看看能不能用，啊这个例子里面好像，啊有一个use multigrade，我可以把它怎么给它换一换，我可以不用毛利柜的，给大家演示一下这个是怎么个传进去的啊。

我把它改成false，我现在把mtv关了，我们来跑一跑，解同样的一个system，我这次不用毛利贵的，我这次只用裸的cg，你可以看到呀，他这个residual掉落速度明显就变慢了。

这已经20个a同学生了，还没解完，multique它有一个性质，就是你multiple写的好的话，那你基本上你的residual是迭代一次，少一个数量级，迭代一次少一个数量级。

所以multigrade它的是这个不管你的这个system有多大，它的迭代次数通常来说是肯定的，比如说呃迭代十次呃，cg 10次，mt grade，那基本上就能够收敛，但是裸的cg啊。

它迭代次数会随着你的，一般来说你可以认为它随着你的信息系统，它的边长是线性的增长，所以cg它的特别是你性能越大，你用毛利贵就收益就越好，你到现在还没有熟练，我就不让他跑，可能跑到下课也不能收敛。

ok大家可以去看一看这个，这个example其实很简单，这个呃这个example呢用了维度无关编程，大家可以把里面dimension改成two呢，就是一个二维的mt柜，改成三就是三维的mt柜。

当然我刚才提到mtv，它设计非常非常复杂的，这个mtv是一个高度简化版的mt贵，他用的所有地方都是最简单的，方便大家去理解它在干啥，ok那么我们总结一下刚才讲了，刚才讲了这么多呀。

其实都在讲这个欧拉的流体模拟，在通信学中，对于每一个time step呢，其实就两个关键操作对吧，advection和projection，the action呢其实就是移动的流体field的。

说白了就是嗯这边唯一一个关键的事情，就是说你要用一个好的advection scheme，使得它的numerical biscoti或者numerical diffusion，要尽量的小对吧。

如果说你的emo with scody很高的话，你这个流体本身不粘，它变得很黏，那你看起来就啊没有那些很刺激的效果，大家在突击圈里面，还是喜欢看雷诺数比较大的流体的模拟，ok然后呢。

如果说你要找一个数值，点型比较小的一个scheme，你可以采用micromi，或者你可以用b f c c，或者你用particle action，后面我们会讲到particle action，第二步呢。

第二关键的东西是projection，projection，就是说你得去保证你的速度场是diversions trade，它得五散，什么是为什么要无伞呢，因为你的素质上五散了以后。

你的流体的体积才得以保持，像我还记得我刚才那个加工吗，如果说我的pressure这个扫不好，如果我只迭代十个一堆人就把他掐了，那我得到的解就不得到，这个压力的解就不是很准确，就会导致速度长是有散步的。

有散度呢，它的这个就不能保持自己的体积，看起来就会呃这个流体会看起来很奇怪，当然了，你要把它迭代到速度上五散，那这个是需要一个很好的零点solver，那通常大家会用mg p c g。

啊对其实我我应该提一下，有一个太极里面有一个sam叫stable fluid。

这个是荒野大神写了一个流体模拟的一个demo。

这个其实是一个交互的demo，大家可以去玩一玩。

这个我没记错的话，应该是用了3d lebron jaction，加上jacob projection，大家如果有兴趣的话呢，可以把它改一改，把它改成microme，然后再改成嗯。

用至少得用contringredient to solve这个pressure吧，当然你如果用stable free和这个呃jacob tation，还有好处，那就是它比较容易实现对吧。

大家可以看到这个还是有一点点黏的，看起来有点黏，然后它的压强也不是完全无损，你如果很猛的去拖一拖，你会看到它这个体积有被压缩或者膨胀的，这个现象，可刚才提到了，其实，欧拉法解这个流体其实就两个东西对吧。

the election和呃和这个projection，然后啊这个vision projection它内部虽然你要做好，但他们怎么组合，这个也很有讲究，这边呢呃介绍一个问题，也就是说你如果原始的速度场。

你做了一次的vision以后，然后你会得到还是这个黑线，你会得到这么个速度长对吧，你得到速度场里面呢，就会有一个有散度的分量和一个无散的分量，这个地方你可以轮子做一个后后侧分解。

然后你这个地方再去做一个projection以后，你就会发现你把有散度的那个速度上，分量给它投影掉了以后，它的动能就会得到损失，它动能得到损失以后呢，那你这个速度长，它就会呃这个地方就会引入一个呃。

numerical viscosity，或者呢你会导致他的vertiity，它的这个涡量它就会给被丢了，被丢了以后，那你就得想个办法把它恢复出来，那这个地方大家可以看张星星师兄的这个文章。

叫做a b u c k是一个挺有意思的文章，他呃通过去重建丢掉了涡量的速度场，给他提高这个流体的这个降低这个流体的数值，连性，另外一个很有意思的paper叫做vtureflection，sober。

他的idea是什么呢，如果说还是刚才那样，我们做一个advection以后再做projection，那么我们就会丢掉这一部分能量，那reflection干啥呢，我去我不去把它嗯。

我不去把它一次把它project掉，我去做一个在advance中间做一次projection，这个projection是能量对于能量没有损失的，我就把这个点给它移到了这儿。

然后后面我们再做一半election，然后最后的这个projection，它能量损失就会小了很多啊，这个这个paper大家可以去尝试实现一下，这个paper非常容易实现，可能五行代码就能实现。

你如果有一个microsoft，再加上一个好的pressure的一个projection，solver的话，可能稍微改改代码，改五行代码就能实现它，而且效果非常非常的惊人啊。

这个是sopa 2018的一个，我觉得非常好的一个paper，既容易实现，而且效果又非常好，大家作业一样，如果有兴趣的话，可以实现一下，这个我觉得呃收获会很多。

而且可以看看他paper里面的推导还是很有意思，ok那么我们基本上就讲的差不多了，呃如果说大家还想继续深入学习呢，毕竟我觉得这节课其实内容也很多了，这个大家肯定要一段时间消化。

嗯如果说大家觉得有一点有一定难度呢，那其实也没有关系，我在课件里面写了好几本书，然后大家可以看看这些书，然后如果呃觉得大家听了这一讲觉得不懂，也不用特别觉得frustrated或者什么的。

不用特别灰心丧气啊，呃一开始总是会有一点难度的，我一开始学的时候和大家其实也一样，但是我只是呃我在这一讲里面，其实更多的是一个导游的作用，更多是一个导游作用，给大家推荐几本书吧。

我觉得你真的要一个小时把欧拉法呃，欧拉视角下的这个流体模拟全部讲清楚，还是非常非常有难度的，我真的已经竭尽全力了，呃大家可以去看看那些书，我觉得还是写的很好了，如果大家还想还想继续提高。

然后首先可以去考虑考虑3d的怎么模拟，然后呢怎么把boundary condition给他做的更好一点，我刚才提到这个boundary condition啊。

只是以voxo为精度的bury condition，很多时候大家会，如果说你的这个边界，正好在你的这个boxo内部划了一刀咋办，这个牛也要做去去做这个cos之类的操作。

那么大家可以去看这个chris bi的这个paper，这个是一个很好的paper，然后流体怎么和固体couple怎么去耦合，也是一个很好的话题。

怎么去做这个two face for resimulation，two face free simulation啊，其实最简单的例子就是说呃你如果有一瓶水，然后假设你有一瓶矿泉水对吧，你把它瓶盖打开。

你把它倒着往下，让这个水往下倒，那你这个水往下流出的过程中啊，你还有空气进去，那你如果说要做一个one face fluid simulation的话，那你这个水就会直接的往下流出来，而没有空气进去。

就会非常非常不自然，如果你做two face呢，你就会有空气进去的，这个效果就会非常的非常有趣，大家可以去看reference，11，是一个很很好的一个用stream function。

来去做two face fluciation，那么刚才我讲的其实全部都是，应该是说是可以认为是气体的模拟气体，它有一个特点，大家一般也叫它烟雾模拟，然后如果要模拟水的话，你就得有一个水和气体的界面对吧。

这个界面它是压强的边界条件，压强的利克雷边界条件压强为零，怎么去handle流体界面，怎么去capture，流体界面的变化就要用到level set，level set。

这个最权威的书就是这个level set methods，and dynamic implicit surfaces，作者是senior shell和rog physical，然后大家可以去呃。

这个这个书其实是一个不是特别友好的一本，对于graphic来说不是特别友好，大家可以去看run他早期rptical，他早期有很多用level set的paper发在sram上面。

那些paper是把这个level set的理论在，在重新学里面的一些比较简化版的应用，然后还有vertex method，vertex method也很有意思啊，那就是说我们之前讲的。

全部都是以速度为基础的n s方程，那么以速度为基础，你很容易做到这个动量守恒对吧，你可能会拖到质量守恒，但是很多时候你很难做到涡量守恒，如果你用一个以涡量为基础的这个server的话。

你就比较容易做到模式量守恒，涡量是什么呢，微量就是流体里面打圈的一部分，如果你做了这个微量守恒的话，你就能保持流体和流血打圈的部分，那你就流体效果就会更好一些，这个我们留给下一讲呃。

新兴师兄来给我们做一个呃，专题专题的这个课程，ok如果大家对实现epied oilers free sober有兴趣的，也可以把它来作为homework，he像我之前讲的呃。

这些应这个oiler resolver通常很难呃，比一般来说你如果叫做要做incompressible的话，那一般来说是没有显示格式，因为显示没法做了，incompressible啊。

除非你用vertistic作为你用涡量的方法，但是我们都是用速度来做，那么你如果想用erance loser来作为homework一呢，那就可以进行implies sover内部的比较。

比如说dj kobe versus，这个country regredients或者multique versus，普通的country degredients对吧，ok那么这一讲就到这里。

我来看看大家有没有什么问题，对我应该终于可以开一下我的虎牙的这个窗口，看看弹幕里面有没有什么问题，我之前不敢开，我怕怕他卡，我看看，哦我现在登录上去，看不到之前的弹幕了，哎呀好吧。

那大家还是把问题发到qq群里面吧，发到微信群里面吧，不知道之前有没有问题，之前之前大家应该弹幕里面发了很多问题，但是因为我没有开播放的原因，我没有看到弹幕，所以不好意思，如果大家想有问题的话。

可以在论坛上面或者在微信群里面发出来，ok或者要不我们今天就先到这儿，然后大家有问题呢，可以在课程微信群或者在论坛里面提出，毕竟其实这节课啊确实容量还是挺大的，我们本来是准备讲两讲。

好多人fluid simulation，但是我们这次压缩成一讲，所以给大家吃了个压缩饼干，可能需要消化一段时间嗯，我一开始学的时候也消化，课要也消化了很长时间，所以我也能理解，这个学起来确实没有不能呃。

一下子掌握，我觉得很正常很正常，这个大家可以看一看样例代码，自己去写一写，有的时候你即使你不理解其中数学，你把这个代码拿下来自己改一改，改着改着，你可能也就理解他这个额数学是什么样的。

它背后数学什么样的对吧，你可以跑跑别人的代码，然后看一看，ok那么我们今天就到这儿，然后感谢大家来我们的课程，希望大家呃做出更加惊艳的作业，我觉得这个alern foresignation。

能做出很多很棒的效果啊，特别是烟雾啊，流体啊，烟雾啊，液体啊，还是很容易做出很炫酷的效果，ok那么就到这了，拜拜我们下周再见啊。

GAMES201：高级物理引擎实战指南2020 - P5：Lecture 5 多体问题与涡方法 - GAMES-Webinar - BV1ZK411H7Hc

那么欢迎各位同学啊，我们今天请来了一位重量级嘉宾，就是我们的张欣欣师兄，然后张毅师兄是呃u b c的计算机图形学博士，然后他给自己的描述是热爱物理与计算，数学，解决其中有挑战性的问题。

致力于在计算机上再现大自然的美，真的是这个，我觉得在座的各位应该有很多同学都很喜欢，希望在计算机上再去模拟一些，大自然里面的自然现象，我觉得呃在计算机上再现大自然的美。

应该还是一个非常非常非常有趣的事情，这也是为什么我们要开这门课的，初中说到这个在线呃普及物理型的各种工程啊，这个其实我之前和呃星星师兄，有一次聊到这个事情，然后我当时有一个困惑。

我就觉得如果我们开了这么个课，会不会有同学觉得我们在拔苗助长，或者呃会会不会有同学觉得超纲了，然后进行实用，就反问我一句，你自己不就是这么拔苗助长长出来的吗，我想好像确实也对啊，只要大家有兴趣。

这个课程稍微难一点啊，或者呃遇到一些难度啊，其实在兴趣的驱使下，我觉得这些都不是问题，当然我们当然会尽我们最大的可能性，把这个课程以一个通俗易懂的方式给他讲出来，ok那么对另外一个事情就是说。

我们以后还是希望大家能多看看直播，因为我们确实有同学反映，我们之前直播有卡顿啊，或者呃背景有噪声的问题，但是我相信现在都已经解决了，因为我换了一台比较猛的g比较猛的电脑。

然后我还专门去amazon买了一个，应该还算比较靠谱的这个microphone，我希望希望以后就不会有这个卡顿，或者呃背景噪音的问题了，而且大家如果看直播的话，有很多问题就可以在线提出来。

然后我们收到这个反馈啊，其实也给我们讲课的过程，其实也是一个非常非常大的帮助，ok那么我们就把剩下内容交给星星师兄了，星星今天会给我们讲讲补充方程，与其快速解法，all right，so，那我们开始吧。

好的好，谢谢这个胡世纪的介绍哈，那个呵呵呵心流水，可见这么多天的授课，这个叫什么口才方面已经长进很多啊，那个对，我希望今天的课尽量能够做到通俗易懂，这四个字，这个但我也相信以我看了上一次的课件哈。

我觉得能够留存到今天的同学，应该应该挺不容易了，所以所以哈那我呵呵，我也不多说了哈，我们现在就开始吧，今天我会讲这个有关播送方程和快速解的，他的一个一个快速解的中呃一些一些事情。

那么这个我们从快速方法来看。

这个泊松方程的时候，我是基于这个泊松方程和它的基础解，这样一个角度来看的，也就是说我们现在讨论的是，把之前我们看到的解法，都是基于它是一个偏微分方程的角度去看的，那么今天呢我们考虑的是一个开空间的。

泊松问题，并且把它通过这样一种积分积分，这样一个基础解积分的形式，来看待这个活动问题，那么其实像这样的一个播送问题呢，很多同学他在你们应该在这个高中的时候，就已经接触过了。

比如说我们十分熟悉的这个万有引力，那么万有引力呢，我们通常就是说我们的市场满足一个，关于密度的泊松方程。

然后呢我们那不同的任意两点之间的受力呢，就是这个市场的梯度，那么在这样一个情况下呢，我们如果把这样一个梯度运算代入到，我们的这个基础解的这样一个求和系统中，我们就会把就会得到我们高中的时候。

很熟悉的那一个呃万有引力的公式，我们的利与距离的平方成反比。

当然这里因为这还有一个。

这个是这个距离的这个方向，所以他一除是一个距离的平方，我们力与距离的平方成反比方向，沿着这个粒子与粒子之间的，这样一个相互的方向，并且前面还乘一个系数对吧。

它与质量成正比，那那么这个呃。

当我们看到这个呃求和表达式的时候，其实我们就能够注意到哈，如果我在一个万有引力系统中有n个例子。

以及我有我要去对m个点求这样一个力的话。

那么我的一个直接求和的办法呢，需要n平方的时间，或者说在这里是n乘以m的这样一个时间，那么这个时间呢它的这个skin是很差的，所以我们需要快速的方法，来降低它的这样一个skin。

那么呃这个时候呢，我们就可以看看一下这张大图，我之后会在之后的这个slides里进行一个详解，那这样一类的方法呢就叫做快速多机展开，这个呃胡适弟应该也在之前的这个呃斯莱呃，之前的课程中有讲。

这个是被誉为20世纪，21世纪十大重要的算法之一，的快速波及展开，那么快速突击展开，它究竟利用了怎样的一种思想，他究竟采用了怎样的一个分析的方法，使得最后人们能够在线性的时间里去计算。

这一个n平方的求和问题。

我们现在就来展开一步步分解看一下，首先我们考虑一个二维的，为了为了简化，因为我尽量希望使这个概念可以清晰的传达，而不是它呃其中具体的一些细节，那么我们现在先来看一个二维的播送方程，和他的一个呃。

这个和它的一个基础解，当然这里为了方便我采用负数来表达，我们的这个坐标，并且呢呃我们知道这样一个球和表达式，它的实部就是呃我们要求的这个呃波动方程，它的解，ok而这个log of距离呢就是二维情况下。

我们的格林函数，那么此时我们考虑一个考虑一个呃，考虑一个圆，以及它在远处这个z点处嗯，引发的这样一个一个一个势能，那么这个时候呢，如果我们在pc on的这个地方去做一个泰勒展开。

那么实际上我们就会发现其实在这一点，在这个z点，他的这个是其实是可以看，至于qi这个这个电子的带电量在零点，在这个零点处，在这个零点处对它的一个势能加上一个高阶项，加上一个高阶项。

这个是我们对于对于对于原来的这个函数，在零点做泰勒展开以后直接得到的这个关系，那么这个时候我们想象，如果我在靠近零点的地方，有很多很多的带电的电荷，它对于这一点的作用。

如果我们都去对零点做一个泰勒展开的话，那么我们实际上就会得到什么呢，我们的fez是等于所有这些东西的，它的这个电量的核对吧，它的视作用的一个和，那么我们在这里对零点做了泰勒展开呢。

我们就发现它其实可以拆成两个部分，第一个是它就好像是一个来自于零点的，一个带了更大电荷的电子，对于这一点的这样一个视作用，以及第二个我们的一些高阶项，我们的一个高阶项是来自于对它的这样一个形，状。

所以这个时候我们看一下，我们把这个公式抽象一下，我们就会发现我们会有一个大q，这个大q呢，它其实相当于是说在零点附近，所有的电子的电量的总和，以及我们的高阶项q k，他呢是在这个零点附近的。

做完泰勒展开了以后，我们得到的这样一些高阶局，像而我们关于这点的这个事，就所有的这些呃零点附近的电子，关于这一点是，它其实是可以写成q乘以我们的格林函数，以及加上这样一些高阶项。

除以z k的这样一个这样一个累加，那么这个后面的这个高阶项呢，我们会发现它是几何收敛的，那么实际上在这个时候呢，我们已经可以得到一个n乘以logn的算法了，这个算法呢又被称为这个吹口虚构的。

他通过这样来实现，首先我们假设自己有很多的电子，在一个空间里头，然后我们对它做一个网格划分，然后呢在每一个格子的中心点，在每一个格子的中心点，我们去近似的计算这个q和q k也不是近似。

就我们去计算q和q k使得什么呢，使得离这个格子远处的这样一个市场的，是可以被q和q k，近似的，那么这个时候你可以想象我们不断的往上，不断的逐层往上抽象，我们最后就得到了一个特别特别粗的一个网格。

它里面存着这些q和q k，其实是代表着所有属于这个格子的这样一个电，量，所所代表的这样一个泰勒展开的这样一个形式，对吧，那么这个时候就有意思了，如果我们想要求这个空间中。

任意一点的这样一个电场式或者重力势，那么实际上我们可以怎么来做呢，首先我们通过树状结构的便利，走到最后的这个叶节点，然后我我们可以，我们认为所有离这些叶节点最近的格子，我们通过直接求和的那个公式来计算。

然后当我们逐渐的这个距离变远的时候，因为原厂它基本上是比较平滑的，那么这个时候呢，我们就可以用那些记录在这个格子中心点处的，这些泰勒展开，去近似他的这样一个原厂式的影响，然后直到最后我们通过最粗的。

最粗糙的那一个解析度上的q和q倍，去近似他的这样一个原厂式的影响，这我们就得到了一个对于这个格子十分精准的，这样一个原厂式的计算，而这个过程呢，因为我们要对每一个格子都去做这样一个，落根的计算。

所以最后我们总的计算量是n乘一老根的，那么我们这边我觉得是有一个问题，就是说包括我也有这个问题啊，就是说这种类似于像八叉树的这种近似的方法，有的地方叫它吹code，有个地方叫它这个burnhot。

然后还有很多很多其他的名字，他们都是一个东西吗，是一个东西全都是一个东西，事实上这个缺口就是bernhart发明的，或者发现了，应该说是对对对，所以这些说这都是一个东西，对没有问错。

ok那么我们更进一步的，我们在知道了这个n log的缺口的了以后，我们怎么样去得到这样一个线性时间的算法呢，这个时候我们就再往前考虑一步，而这个问题就是我们能得到欧文时间，算法的一个关键。

此时我们假设我们已经在一个z1 ，这个中心处的电子云，我们已经统计了出它的这样一个q和q k对吧，它可以精确的估计，在原厂z处的这样一个电子式，那么这个时候我们要问一个问题。

我如何通过这个电子云的q和q k直接估计，估计出另外一个z2 处的q和q k，使得在这个z点处的，是可以被z2 处的这样一个展开，所精确的表达呢，这个时候我们其实涉及到的就是什么呢。

就是把一个multiple，把一个这个multiple呃，所以这个又叫做一个m2 m的变化，那这个时候其实方法也很简单，我们既然有了z一处的这样一个电子式的，这样一个计算公式。

我们对它重新在z2 处进行一个泰勒展开，并且重写这个公式，我们就得到了一个新的这样一个形态，你会发现，其实在这以z2 为中心的这样一个，这样一个电子云，它对最初的影响可以写成这样的一个公式。

就是q乘以log z减z2 ，加上一个suv加上后面的高阶项，我们这个时候recall一下，我们刚才得到的这个multiple展开的公式，我们发现其实他们俩是有一个类似的形式的，首先我们都有一个大的q。

代表的是所有电子电量的和对，然后呢我们都有一个q k，当然在这里它变成了一个pk对吧，代表的是一些高阶的余项，而他们的收敛性呢也是都是这样一个几何收敛，那么这个时候我们就说dk到底是什么。

其实bk呢我们发现它是对于q k的一个泛化，它是对q k的一个泛化，什么叫泛化呢，首先我们来看bk的第一项，它是大q乘以z一到z2 的k次方除以k，这个东西是什么呢，我们recall一下。

在之前这个我们在看算这个，我们recall一下，我们之前在算q k的时候，他的第一项是什么，他的第一项是关于这个qj的一个和，ok，那么在这里呢它就变成了这个他的第一项，也是关于qg的一个和对吧。

所以其实这两项是相等的，ok那么以及呢它还有自己的一个高阶项，这个高阶项是通过这里的展开了以后，直接这个推导计算出来的，ok，那么在这个时候我们发现它其实这个高阶项呢，是关于原来的那些qk的。

那些高阶turn的这样一个计算的表达式，那么它实际上形象的来看，这相当于是个什么呢，我们的这个我们的这个电子的电源，它其实就是一种十分特殊的model，他们带了电量，q which当它是一个电子的时候。

他就是他自己带的电量，他的所有的qk项都是零，而这个时候呢，我们直接套用bk的这个计算公式，我们其实是可以复现我们的这个multiple expansion，也就是说我们可以复现出从电子累加到这个呃。

累加到一个点的这样的一个multiple expansion，那当我们从n个或者当我们从呃不对，当我们从s个这样的已知的mmultiple，去构造一个新的multiple的时候。

其实每一个已知的multiple，它就又相当于是一个电子了，它就又相当于是一个电子了，然后这个时候呢我们就computer这个bk，通过使用什么呢，通过使用呃，呃加上这些这个rest of turns。

那这个时候我们的我们还是要注意，我们的bk的第一项的这个q，也就是这些小的电子云，他带的这样的一个总电量，所以这样一种操作，我们实际上是把multiple直接给translate。

到了一个新的multiple身上，那么可能刚才讲的有点迷糊，我在这里呢就用一种代码的形式，再把刚才的过程展开一下，首先我们可以定义一个multiple的结构，这个结构呢它包含了两个重要的信息。

第一第一个是它的终点坐标在哪，它的中心在哪，第二个呢是我们有一个，我们有一个qk这样一个turn，当然在这里我们讨论的范围内，它是一个负数对吧，它来表达他的所有的这些高阶项。

那么当我们从source charge来出发的时候呢，它其实是一个特殊的multiple，这个multiple呢它的center就是我们charge的这个坐标，以及呢它的这个q0 。

他的q0 就是他自己带的电量，它其他的这个q k的这些turn呢就是0k，这时候我们可以写出一个通用的multiple to multiple，的这样一个圈子transfer。

那么它呢take的这个参数是一个q list，是一个vector of multiple对吧，然后这个时候我们先去计算它新的终点坐标，比如说终点坐标。

我们可以通过对于这些q list的center的一个加权和，来得到，以及呢我们去计算出它的q0 t跟q0 turn呢，就是所有的q list，它的这个q0 t的和。

然后呢我们对于从一到p这几个高阶展开上，我们呢调用这个，我们调用computer bk的那个公式，去计算它的这个p turn的这样一个值，然后我们把这个result返回，这个时候。

我们就构造了一个从电子到multiple，以及从multiple，multiple的这样一个转变的这样的一个程序。

那当我们知道了如何从multiple，转变到另外一个multiple的时候，接下来就有一个问题变得十分有趣了，这个问题就是，如果我们在一个中中心为c的地方，知道了他的multiple展开。

也就是说这个multiple可以告诉我们，在一个任意一个z点处，它的这个它的这个准确的，这个电场市到底是多少，那么这个时候我们就想问了，我如何能够把它转转换成一个离z距离很小的。

另外一个z一处的这样的一个local口，我们注意这里不是multiple to multiple，transform，multiple to multiple，它的这个transform。

他要求的是我的multiple一和multiple 2，都是离我要求的这个z点是比较远的，这个时候我们是m two m transform，而回到现在的这个问题。

我们希望把一个离z比较远处的这样一个multip，转化到一个离z比较近的地方呢，这个z一为中心的一个local口，这个local local，大家又可以形象的想象它是一个什么呢，它就像是一个差值。

它就像是一个差值，它使得在这一点处的这样一个市市场函数，能够很准确的被这一处的这样一个，多项式给表达出来，所以它就像是一个差值，一个局部领域的差值函数，那么这个时候我们去看一下，当我们做完了。

我们该做的数学了以后，我们其实就得到了这样的一个公式，首先这个z一处的这个这个local它的第一项，它的第一项我们可以看到它是什么呢，它的q乘以log，z一减去c加上这样一个peter的这样一个求和。

它实际上是什么呢，它实际上是电场是在这一处的这么一个势能，这么一个potential值对吧，所以它实际上就是电场是在这一处的值，那ok再加上一个什么呢，我们要求z点数的这个电场式的话。

它实际上就是把电场是在z一处，再加上一个高阶的扰动，其中这个高阶扰动呢由这个式子给出，它是一个bl乘以z减z一的l次方，的这么一个形式，那么由于我们知道这个z减去z1 ，是一个很微小的量。

是小于半径的一个量，而bl呢它的它的半径是这个，它的半径是这个z一减去c的k次方，所以这个这个式子它又是几何收敛，那么这里我们就给出了bl的这么一个表达式，这里我们就给出了bl的一个表达式。

这呢我们就完成了一个，从multiple到loco的这样一个变化，而这个时候我们loco的这样一个结构，它大概是这样的，首先他知道他自己的中心在哪里，然后呢他有一个turn叫做bp。

就是还有这么pig这样的这样一个beer，它来计算，那它的b0 是谁呢，b0 就是它的终点处的电场式的值啊，这个时候我们就又来了一个新的问题，就是如果我们已经知道了，它能够表达所有在它这个半径范围内的。

这些点的这个事的做事的这个值，那么我们怎么给它转变为，一个以c2 为中心的loco，都可以被c2 的这样一个local ball准确的表达，那这个时候，其实我们涉及到的只是一个坐标的变换。

只是一个坐标的变换啊，那么比如说我们就把式子直接写出来，我们其实要求什么呢，是ak乘以z减去z0 的k次方，它可以变到一个bk乘以z的k次方，这样一个形式上，那么会这个可能学完初中的这个数学。

我们就能够做这件事了，那么实际上在这里呢，是有一个这个叫做honor scheme的方法，在计算上帮我们进行一个动态规划啊，动态规划的这样一种计算，可以加速这个turn的这样一个计算的过程。

诶到这我相信有的同学会有一个问题，就是说为什么multiple to multiple，看起来也是一个呃中心的转移，然后locple to local也是中心转移。

为什么local to local就比较简单，而multiple multiple的形式会稍微复杂一点的哦，其实也都差不多是吧，其实也差不多，其实都是我觉得他也没有简单多少，你看这里其实还是很复杂的。

我只是写的简单了，这个形式我觉得是挺复杂的，他的这个公式我觉得其实是挺挺挺烦的，所以我如果不躲都不建议做这个local to local的，我还是觉得用插值的办法会好，构造一个差值就行了啊。

好先这其实它只是形式看起来简单了，我这一个bk对吧，但其实这个bk呢，你会发现它代表的是这么一坨东西对吧，那么在三维的情况下，这个式子会更加的恶心，你看这里这是啥，这bono。

这是binomial coefficians对吧，对对对对，这很恶心的这坨东西，ok嗯那还有一个问题，就是说我们在这说的multiple，和在物理里面的电路机子，它们之间是有对应关系的吗。

是有对应关系的，是exactly的一个对应关系，而且那就是他的说b0 是电单极子，b一是电偶极子，是这个吗，对就是这个机子对，ok所以就是dn就是电二的n次方极子，对你都可以这样想象，其实是怎么样的。

其实在三维的情况下，我们那个我们那个方程是one of them are嘛对吧，是1÷2的那个r i j，然后那个时候我们如果在cartesian coordinates。

在这个笛卡尔坐标系下对他做泰勒展开，然后我们得到了一系列的这个相对吧，这个monopole什么quad pole，八八炮对吧，这些东西就分别对应的是物理上的那样，那样一个概念，ok alexcess嗯。

那么那么在这里其实我们我更喜欢哈，把刚才所讲的这一整个过程再做一个抽象，我们会发现其实multiple expansion，我们实际上在做一件什么事情呢，就是对原来的这个几何，对于这些电子啊。

我我看到我们也把它看作一个几何的话，他是在做一个cosy，他是对他在做一个粗化，然后我们的local expansion呢，它实际上是在做一件什么事情呢，他是在做一个差值。

所以这个时候我们看multiple transform的时候，我们一开始从一些电子开始，然后我们在格子的中心点，去计算出它的qq k和c对吧，然后我们再把格子上的qq k和c。

计算到上一个格子的q q k和c，这个时候呢我们总共的计算量，我们总共会便利过的格子的量，一加上在三维的情况下是1+1/8，加1/64加点点点，所以总共我们会遍历过的格子量是o n。

然后在每一个格子上呢，我们会从事的这样一个计算量，大概是c的p次方对吧，大概是c的t次方的这样一个计算量，所以在这个时候我们往上行的时候，我们的计算复杂度是线性的，那我们再来看下行，下行的时候呢。

我们会经历multiple to locle的transform，以及local to loco transform，那可以当我们在执行multiple to loco，transform的时候。

我们可以注意到，对于每一个格子，对于每一个格子，它所相对应的multiple的数量是有限的，他并不会去累加所有的multiple，他只需要去计算在自己这一层上，和自己相对来说是multiple的东西。

因为来自于原厂，也就是，比如说考虑这个黑色格子，和这些灰色和这些虚线格子，这些虚线格子是这个黑色格子的原厂，这些原厂对它的作用其实通过谁来完成呢，其实通过loco to loco的这样一个变换来完成。

也就是说，当我要计算这个黑色格子的原厂作用力的时候，我实际上是做了两件事，第一件事我从它的复节点，从它的这个复节点去继承loco，并且把附节点的loco给转变到自己的这个，格点中心上来。

然后呢我在我自己相应的这一层上去累加，对于我来说是multiple的这样一个贡献，那么这个时候呢，其实每一个格子它的计算量是常数时间的，而这个常数时间大约也是c的p次方，而这样的格子总共有多少个呢。

一样也是n加上八分之n，加上64分之n点点点，所以呢我们下行时间总共也是欧文的，那么我在这里我再对比一下，这个数代码和快速多极展开，他们两个之间的一个区别，我们可以看到，如果我想对x一和x2 这两个点。

去计算他们的市场值的时候，我们发现对于数代码我们总共要去对于x1 ，我们要累加123456，这些格子对我的一个作用对吧，以及于呢对于这个x2 ，我仍然要去累加123456对我的这样作用。

所以我当我有n个格子要计算我整个长，每个格子呢我实际上要去累加的量，大约是一个logn的计算量，所以它的计算复杂度是n log n的，而回过来我们看快速多极展开，当我要计算x一的所有的贡献量的时候。

其实它总共分两步，第一步是啥，第一步是红框，这个格子也就是它的复节点早已经累加好的，来自于123这三个格子的贡献对吧，这是复节点上就存着的，然后呢我们通过一个l to l的变换。

把副节点上的贡献直接先变换给我了，直接先拿过来了对吧，然后我再去累加来自于这个456，这三个节点的这样一个贡献，那么这个时候对于我来说，我每一次的计算时间，计算复杂度都是一个常数的，所以这也是，这。

也就是为什么fmm，它能成为一个线性时间复杂度算法的原因，那么实际上在我看来哈，前面的这个这些local people do to loco的，这个变化是十分复杂的，我们可以考虑一个更简单的情形。

如果我在附节点的时候，我不去计算local和我取而代之的是什么呢，我在复节点的这几个节点上，这几个box的这个节点上去，计算出来自于123这三个东西的一个贡献，然后我干嘛呢，我在计算x的时候。

我通过一个双线性插值，或者在三维的时候，一个三线性差值先去插值出这个贡献，我不就等于已经拿到了，来自于123的一个贡献吗，因为足够平滑，所以我认为差值也是可能的，ok我再说一遍。

我们在这个节点处去用multiple的方法，去计算1123处的这样一个贡献，然后在里面呢我先插值，插完值了以后，我再去计算456的贡献，并且我把这些东西都计算在哪呢，都放在我这个格子的格点上。

我认为这也是一种方法，这可以避免这个locpto local的这样一个，比较复杂的计算量，比较复杂的这样一种计算，其实我觉得那个代码写出来也很恶心，可能有很多bug都不好修，因为这个公式确实太恶心了。

可能换一个人都修不了，修不了bug，对啊，我昨天还打算自己看看，一天之内能不能实现一遍，然后写着写着发现了这些公式实在太多了，这个一天肯定写不出来，还是很复杂的，所以大家一定要知道，这是一个。

连胡渊明都不可能在一天内写完的东西，我以前玩的东西太多了。

对哦好，所以我们再来回顾一下这张全图，这张全图其实给了我们一个很好的这样一个后，羿面前，他表述了说，首先我们从这个source party开始对吧，我们先做p two m变换。

然后呢我们再做m two m变换，然后我们再做m to l变化，然后我们再做l to l变换，最后我们做l two p变换，然后呢在particle周围呢，我们再去累加来自于其他粒子的这些贡献。

这就是整个这个fast multiple它的一个一个一个走法，先是top，先是button up，然后再top down，所以其实我们会发现这个东西的计算。

它似乎和multi grade method它也有异曲同工之妙，也存在这样一种空间的这样一种呃，button up和top down的这样一种结构，那么在三维的情况下，在这儿我要指出一下。

就是我们的泰勒展开可以以两种形式来完成，一种呢是在ktition space，在这个笛卡尔坐标系下进行展开，另外一种呢就是在这个球鞋坐标系下，进行这个展开，那么在球鞋坐标系下进行展开的时候呢。

我们就得到了原始的fast multiple method，那这个大家如果想看更多的这个细节呢，我推荐大家去看这个a short horse on，fast multiple methods。

因为他是这个green guard写的，green guard呢，他就是fmm这个算法的一个发现者，所以这个这个课程我认为是呃，这个cos我觉得是写的比较好。

有同学问，是不是任意的是函数都可以用多级展开。

这个就有意思了，这个有关于什么样的势函数。

就是说有关于什么样的和才能够被施展开，在这一在这个short course里面，它是有有讨论的，他有一个引理的证明，就首先你这个和是要满足一个条件，你才能进行市展开的。

那么这种核呢我们叫做degenerate kerne，嗯哼ok它要符合这样一个条件，那么通常的物理上的来来，源于物理应用的这些和，我们可以认为他都是可以被这个市展开的。

包括高斯和包括这个呃热热力学的这个和，包括我们的这个呃格林函数啊，这个万有引力的这个和它都是可以被展开的，那么至于任意的和能不能够被展开呢，大家就要去套那个引力验证一下啊，这个不是任意的，ok政治是。

看来这个课程这个show cos确实非常给力啊。

我也应该去学习一下。

那实际上这东西呢它这个fast multiple method，它是有很多这个呃有趣的应用的啊。

一个比较重要的应用，比如说我们可以计算gravitational force。

以及关于重力的场的一些信息的变化。

那个这个除了可以帮我们在cosmology这个方面，预测星体的运动以外呢。

我们还可以就还有人用它来搜寻这个dark matter。

那事实上这dark matter这个名字很酷炫啊，暗物质或者暗能量啊。

那么嗯如果我们搜这个fast multiple method的话。

我们其实是会搜到有关于dark matter的，这样一些论文的。

那么它实际上它是要解决一个什么问题呢，当我们观测到一些带有星体的质量。

或者我们对一块区域有这么一个质量的估计。

以及我们知道它在怎么运行的时候，这个时候我们就可以借助这样一种。

仿真分析的办法，去尝试去搜寻出我们计算出来的这个轨迹。

是否真的符合我们观察到的现象，如果这其中差异很大的话。

这说明什么，这说明就存在我们没有看到的质量。

那这不就是所谓的暗物质吗，啊所以如果对这块有兴趣的同学。

也可以去看一下有关于这方面的论文哈，那么另外一个它比较有巨大作用的领域，是什么呢，是静电场，electrical static，electrical static，它和我们的这个万有引力一样。

它具有完全一样的这样一个泊松方程，当然正负号不一样啊，这里我要强调一下，正负号不一样，但它的形式是完全一样的，这样一个泊松方程，它除了可以在这个静电物体，物体的表面就给定一个电子的分布。

去计算出静电场的这个静电场式的，这样一个方向以外呢，它还可以干什么呢，事实上在分子和分子之间，或者在分子尺度，静电作用力是它的主要的一个作用力，那么这个东西呢，就会被应用在分子动力学分析中。

那事实上在一个分子动力学软件中，每一个时间去计算这样一个分子静电场力，是呃十分核心的一个任务，而对于分子动力学的分析呢，可以有助于人们做这个肿瘤方面的研究，以及做这个靶向药的设计。

以及像现在呃像呃新冠疫苗的时候，人们就用这样的方法来分析这样一个新冠疫苗，它最好的这个a c r受体的结合点的位置啊，这其实可能很多呃，大家可能很多都不知道，原来这个之中还有这样一个泊松方程。

在再play一个important ro。

长见识了哈哈，那么以及以至于我们有这个hao方程，hamhz方程呢，它就通过我们对这个hamhz方程的，这样一个求解的，我们就可以去求解这个生产，比如说在这个声音领域。

当我们在某一个位置放了一个这个扬声器，我们想计算这一整个范围内的这个整个声音的，一些声压的分布等等等，这样一些情况呢，我们就可以通过求解harmless equation，来来解决这个问题。

那么以及于呢它可以在这个electro magnetic，就是电磁立场这样一个方向呢，呃产生一个应用，比如说当我们在设计某种健体的时候，那么我们就可以通过这个电磁力的这样一个分。

析去找出它可能被雷达探究的半径，可能被雷达探测到的半径，以及我们去反向优化它的这样一个形体设计，使得它尽可能的能够，呃这个达到一个更好的隐蔽性。

那么这个就是通过electro magnetic的计算得到的，而这两个东西它都来自于这样一个ham，ho或equation的这样一个求解。

那么以及以及呢，我们有一个十分重要的应用是什么呢，就是这个边界元方法，在边境源方法中我们求解一个什么，求解一个拉普拉斯方程，这个拉普拉斯方程呢，它其实是告诉我了我们这个整个流程，整个厂里头是什么。

是有事无缘的对吧，有事无缘，那么这个时候就有意思了，其实这样一个有事无缘的厂子，它的所有的属性就被边界出的值，所唯一的决定了啊，所唯一的决定了，比如说我们如果在边界边界处。

这个比如说在这里这个伽马一这条边界上，我们知道了他的def dn，以及在伽马二这条边界上，我们知道了他的b的这个值的话，那剩下的整个厂里的所有的数值，就都是确定的了，这个数值的确定呢。

可以通过一个解析的表达式来得到，这个表达式呢就是边界积分方程啊，我这里已经把它写成离散形式了，我这儿已经把它写成离散形式了，所以实际上我们是什么呢，我们把一个偏微分方程，带边界条件的偏微分方程啊。

转化成了一个积分方程，在这里我们把它转化成了一个积分方程，这就是边界元方法，那么这个积分方程它有一个特点，就是说它的矩阵是稠密的，因为每一个点它都来自于其他所有点的作用，所以它的这个方程是矩阵是稠密的。

但是usually它的肯定是number是很好的，比如说如果我们用一个co subspace的这样一个，迭代法去求解这样一个积分方程的时候，通常我们在常数迭代时间里就能够找到他的解。

所以他的肯定是number，通常是很好的，但他举证很稠密，所以这里我就列下对比，在二维的情况下，在二维的情况下，假如说我们fuldomain是n平方，负的domain是n平方，这个时候。

如果我们用multigree的方法来求解偏微分方程，由于multi的方法是在常数时间收敛的，那么我们得到的是一个，总共是大概一个n平方的一个计算量对吧，而如果我们用来if。

或者用strike forward，的一个边界元方法来计算的话，那么我们要知道边界圆在一个n平方的范围内，因为它只有边界嘛，所以大概有n个element对吧。

那比如说我们用by cg step这样一个方法来迭代，它可以在常数时间收敛，但是每一次迭代，那它要用n平方的时间来计算，这个稠密矩阵相乘的这个时间，所以他最后的总的计算量呢是和multi gi的差不多。

打平手的，就如果他买衣服的来实现的话，他是打平手的，那在三维的时候情况就糟糕了，比如说我们的fdm是n的三次方，那multique的方法呢。

我们知道它是o n3 的时间收敛，那么当边界元呢，这个时候他有n平方的elements了。

而计算这个从密矩阵的代价呢，就是n的四次方的computation了，所以即使它能在常数时间收敛，它总共是n的四次方的computation，它就和这个multigree的方法。

他就永远都无法战胜multipi的方法。

那么反过来，如果在计算这样一个积分求和的这样一个时候，ok计算这个积分求和的时候。

我们用之前所推导的快速多极展开的方法，用快速求和的方法来替换掉，这个矩阵向量相乘的过程，简单的说原来它是一个矩阵乘以向量，这个时候呢我们用matrix free，用matrix free的这样一种形式。

然后用fmm去计算它的这个结果的话，那么在2d我们就是o n的，我们打败了多重网格方，然后在3d的情况下呢，我们就是o n平方的o n o n平方的，我们也打败了这个多重网格法。

那么我们要知道它有一个好处，就是边界圆这个方法它是sumalthletic的。

它是一个半解析的这样一个方法，在我们对于圆的这个划分，精细到一定程度了以后，他给出的答案几乎就是解析解，所以这个是他十分有趣的一个性质哈，他给出的答案几乎就是解析解。

并且很呃不用很细致的一个这个边界源的划分，就已经可以得到很好的答案了。

所以这个是他很有趣的一个性质，它是3m ltc的。

那么比如说用这样的一些编辑的方法，我们可以来求解哪些问题呢，那ok比如说在这里我就展示了说，有关于波动的这样一个这样一个呃方程，那么这个是像像我现在所展示的这样一个，simulation domain。

如果我们用呃n s方程的方法来求解，我需要把整个抖妹抖体积化，这里面就需要很多很多的计算员，这里面就需要很多很多的计算员，而我如果用边界圆的方法在上面求这个呃，式方程，或者说叫做伯努力方程啊。

我可能得到的物理，没有n s方程所表达的那么复杂，但是要知道这个呃国动力方程，已经是对n s方程很好的一个近似了，因为他只是假设你的流畅无年而已，只是假设这个流畅无年。

而且流畅的这个速度是一个是一个市场，是一个市场，那么在这个情况下呢，我们是可以求解出复杂的这些波浪运动的，这样一个形态的，那么对于很多应用来讲，由于波浪的这样一个，特别是远处的波浪的这个这样一个形态。

已经是很核心的一个研究任务了，所以在这种情况下，使用边界元的方法求解，很显然是更明智的，那么以至于在这里所展示的这个案例给我们，我们可以去计算升压，通过求解之前所说的求解这个emo equation。

这样一个边界圆方法，我们就可以说，在我们已知一个物体的外形的时候，我们就可以在它周围，任意范围内去计算这样一个升压，那计算升压有什么好处呢，我们就可以知道它在飞行，在这样一个以一定速度飞行的过程中。

所产生的噪音，或者所产生的可被探测的这样一个声音，它的影响的范围是多少啊。

我们可以去进行这样一个升压的计算，那么以基于什么呢，在图形学或者在很多这个有关领域吧，我们还很关心的就是湍流，也就是涡流的这样一种计算，那么涡流实际上呢，我们来自于一个矢量的泊松方程。

我们首先通过求解一个矢量的泊松方程，来得到我们的一个呃string function流函数，然后我们再对这个流函数打这个客运算，来得到我们的速度差，那么它反过来呢，又能够写成这样一个求和的问题。

那么这个时候呢，我发这个fast summation，又可以用来求解这样一个涡流。

那么求解涡流这个问题。

他其实特别有意思啊，比如说我在这里展示的这个案例，我在这儿展示的这个案例。

所所展示的这个现象叫做vtex ring。

他是讲这个两个vtex ring，一旦产生运动了以后。

他会做这样一种你中有我，我中有你的这样一种交叠上升，ok那么这样一个这样一个现象。

如果我们想通过n s方程去捕捉。

其实他是十分痛苦的，比如说我在这里需要用很多的这个自由度。

我用了一个512的平方的这样一个网格，然后呢我要用一个advanced呃。

比较高级的这样一个时间积分的方案。

高阶的这样一个时间积分的方案，并且这个server呢，他可能要很多很多行的c加加代码。

即使用太极来写，它可能有很多行。

它才能差不多的模拟出这样一个leap。

而最后我们发现他的这个box还是，很快的就融合了。

很快的就融合了，那么现在我们看一下，如果直接用这个握方法来写的这样一个东西，我们我们就可以捕捉出这样一个leap，他呢我的时间积分方法，我可以as simple as s一个前向欧拉。

并且我在整个计算中只有四个自由度，就是四个我点那四个窝点，这四个自由度ok也不，当然也不止四个4x3，他们的坐标对吧，这样这么点自由度，反正然后呢我只要很少行，我记得应该是数应该是百来行的c加加代码。

那么在太极里头呢，它的代码的行数就更少了，它大概只需要几十行呃，多少行，七八十行嗯，居然也要七八十，你是算上visualization了是吧，嗯对，就是其实里面好像我记得，好像初始化这两个窝的位置啊。

还有什么的，就花了两个10号，哈哈行吧，就就刨掉这些以外，大概只需要40行左右的代码部分，大概40吧。

对，我们就能够捕捉出这样一个复杂的，这样一个运动现象，而这个要知道我们是通过n s方程的求解。

是要花很大的力气才能够做到的。

所以到现在我差不多总结一下，我们总共有这么几种submission是要求解的，比如说第一种关于电场力和万有引力，那么就是市场的梯度，市场在x点的梯度，那么它的求和是关于平方，关于半径的平方成反比。

质量成正比对吧，方向是沿着这个距离的方向啊，这里消掉一个消掉可以消掉一个r啊，那么以及在这个边界园中呢，我们要对一个不同一点的kernel做计算了，这个kernel是说我需要在每一个j。

每一个每一个便利的j的位置去点击，去拿我这个方向的这个梯度，去点击我当时的那个法线，点击j a这个地方的法线，那么这个会使得这个kernel变得很不一样啊，这会使这个ko变得很不一样。

所以呃我们要考虑一下这个东西，我们怎么来求对吧，以及呢我们有一个bo smart，bo smart，当然他很像我们的这个万有引力的，这样一个公式了，它只不过这里用的是一个叉乘的形式。

ok大概这三种submission。

那么这个时候我就有一个问题了，我就说如果我给大家一个routine，这个routine叫做computer graduate，如果我给定一些电子元或者重力源s。

这是n个particle with its mass，对吧，这是n个particle with its mass，然然后呢，我来computer grab fee，我就能够告诉你。

所有的这些po位置出的这个gravy，也就是他们的受力，这个f的话，我如果给了你这样一个routine，你把它当做一个black box。

我们怎么去求剩下的这两个，我们怎么去求剩下的这两个求和。

那么这里我展示一下，我们怎么求这个边界积分中的，这个nj的这样一个求和，其实我们可以把nj有关于nj的那个东西啊，我们把它拆成三个特啊。

我们把它拆成三个ten，我们发现其实是ngx乘以d，这个one nor d x加上ng的y分量乘以d。

他d他的这个y分量，那么这个时候我们就我就发现这件事情，他就很有意思，他可以这样来做，如果从原来的这些s，从原来的这些s我们去构造出三个不同的s，构造出三个不同的这个圆，那么其中呢这些圆这样来构造。

我把每一个si自己的原有的这个带电量去乘以，他当时的当初的这一个法线的x分量，或者y i分量或者z分量，我去构造出三个新的s，并且呢，我对这三个不同的s都去做一个computer grade feed。

这样一个操作，ok然后我分别把这个f一的x分量，f2 的y分量以及f3 的z分量，我把这三个东西加起来。

实际上我就得到了，对于这个式子的一个求和的值。

我就得到了对于这个式子的一个求和的值，那么中间具体的这个呃，中间具体的这个过程，我希望同学们可以自己拿个笔写一写啊，你就很容易能够想通这件事情，这是很显然的，那么至于bo smart怎么来做呢。

我就留给大家作为一个家庭作业了。

好吧，在给定grade feat，计算grad feat的这样一个基础上，怎么去计算这个bo smart啊，注意我们是对于j去便利的啊。

注意我们是对于j去便利的，对大家如果做了这个bo smart这个作业的话，也可以在论坛上面和其他同学分享分享，这应该是一个祭祀，主要是一个数学公式的推导对吧，对，主要是个转数学公式推导或者转化的过程。

实现起来还是很难的，我觉得实现起来是不用实现，就是推一下就行，拿支笔一写就行，那么剩下的我呢讲一下其他的两种，这个快速求和方法，一种呢就是p3 m。

一种呢就是这个呃和无关的kernel independent呃，快速求和，那么首先我先讲一下p3 m，其实pcm他的观点就比较的怎么说呢，比较的这个嗯比较的印q一体。

我们是通过这个去混合了偏微分方程的视角，和samson的这样一个视角，来得到p3 m这个方法的，那么我们注意到每一个这个submission，其实它都是对应着一个这个偏微分方程的，每一个扫描线。

它都都对应了一个偏微分方程。

那么这个时候我们就可以，就可以先想这样一件事情，如果我有一些例子，我先把这些例子做一个p to g的操作对吧，我先把例子speed到这个网格上去，然后我在这个网格上呢去计算一个什么。

去计算一个泊松方程的解，去计算这个泊松方程的解。

得到我的流函数或者得到我的势函数对吧，得到流函数或者是函数，然后我对它施加这个梯度或者克，我实际上这个时候就得到了什么呢，我就得到了一个平滑的市场，我就得到了一个平滑的场对吧。

这个厂注意了是包含所有粒子的贡献了，这个要注意，是包含所有粒子贡献的一个平滑的这样一个厂，它包含了，但是呢这个场这个我们得到的场，他的进场的这个值它是没有细节的对吧，它是没有一个细节的。

所以接下来我们要考虑的是，如何把这些细节给加回来，那么加加回来细节的方法就很有意思了事，事实上我们可以先把这个平滑的场插值出来，然后呢再对于周围去累加一个呃，去累加它的一个高频。

高频的这样一个像所带来的影响。

我们就可以加回这个湍流场的这样一个细节，那么在我把这个算法给奥特曼之前呢，我先呃向大家展示一下，那这样一种做法，它首先它的计算时间它确实达到了线性，比如说在这里这条蓝线是o n哈。

然后这条绿线它是o n平方，所以计算时间呢我们确实达到了这样一个线性，那么我们的计算精度呢，我们也可以看见，随着我的correction range的这样一个增加，它能够到达这个小于1%的这么一个状况。

在这里我记得呃，我在这应该是小于0。01%了，所以是达到十的-3次方这样一个精度的。

那么呃，那么我现在把这个p p3 m给奥特曼一下，它其实是这样，首先我们有一堆例子对吧，我们去做一个particle to bread的这样一个操作，然后此时我们得到的是一个平滑层，我们得到一个平滑层。

然后呢我们去计算它的高频场，计算高频场，我们把每个例子减掉，他从平滑城里差值出来的那个值，我们就得到了它的什么高频分量对吧，高频分量，这个时候呢我们拿我们的这个泊松方程。

我们拿我们的这个psp去求解这个平滑的场，并且得到我们的流函数或者我们的势函数，然后我们对流函数或者势函数求curl，求grading，我们就得到了我们的平滑的这样一个速度分量，速度分量。

那么这个时候我们最后在粒子上的这个，湍流的速度分量，它由两部分组成，首先我们去插值平滑的速度分量，然后呢我们再对它附近区域的这个高频分量，这个高频分量去做一个什么，做一个快速求和，也不是说快速求和。

是说叫做直接求和，那么我们就得到了这样一个湍流的，这样的一个速度差，ok所以那这样的一种计算呢，由于我们每次我们所有的计，在整个过程中的计算都是线性的，所以最后我们得到的是一个什么呢。

也是一个线性时间呢求科算法，并且它比这个快速多极展开要简单得多，同学们其实甚至用有用呃，课程至上一节为止，已经开发的这些呃工具和能力，应该就足以写出这样一个呃p3 m算法，这里面这个第三步是用mt贵的。

应该就可以吧，这就用猫体贵的，就用猫体贵，但是你边界条件是不是有一个。

还有个8号出院那边诶，不用你看我这里p3 不用边界条。

就是说边界条件可以认为是零，你看我这个图我画的很清楚，所有的这些例子，它是split到很小的一个格子上去的，很小的一块区域上去的，我就边界足够远吗，然后这个边界我就是让他足够远，就是让他足够远。

因为什么这样可以呢，因为我最后有那个高频分量的一个直接累加，对吧，所以这样是不会损失进度的。

这样是不会损失进度的，对，那么还有一种呢，就是这个和无关的快速求和方法，那和无关的快速求和方法，实际上他的视角是这样的，当我有个这个电子的时候，我去构造两个圈，把它给包住，我去构造两个圈，把它包住。

这两个圈是这样的，首先我来自于呃，我我要在pj这些点上，也就是说在外圈上去求出来自于qi的，这样一个直接求和的贡献，这个是我已知的要达到的数值，与之要达到的数值，然后呢。

我假设这些数值是由来自于内圈上的，一些未知的电子的电量来给出的，这个时候我要求说什么呢，我内圈的电子的电量要是多少，我才能够在外圈上给出这样一个四分布，所以所以而这个k k j是已知的。

这个pj是已知的，然后我们要去求这个西格玛k，所以这个时候其实我们是一个什么，是一个矩阵求逆的过程，locally的一个矩阵求逆的过程，而一旦我们这个内圈和外圈这个形状是确定的，这样一个情况下。

我们的这个这个是这个kk j，就是一个已经确定的值了对吧，所以每个圈上有多少个多少个电荷，这个是随意，这个是多多益善的，就是这个电荷越多，对这个电荷越多，你可能越准确嘛，对吧对。

所以这中间差不多就是有几百个呃，有时候100个吧就够了，ok，啊所以这里我再做一下总结，今天呢我们讲了几种这个快速求和方法，有快速舵机展开和p3 m，然后我们呢我后来讲了这个。

所以能够被快速求和方法来求解的方程，它们包括了泊松方程，拉普拉斯方程，pho方程，还有边界元方法对吧，然后以至于呢这些快速求和方法，它可能产生的应用包括了在静电场方面，它有这个分子动力学。

它可以作用于肿瘤呃，药物设计等等，以及我们的electro magnetics，它可以有这个船舶的这样一个设计对吧，然后我们的这个声声学方面呢，我们可以做这个urban planning。

比如说你现在要建一个高架，这个高架呢会经过几个居民区，那么这个时候呢，实际上可以通过一个升压的分布，去在高架上面呃，适当的放一些隔离带，或者叫做隔离的那个隔离板，来起到一个最好的去弱化这样一个呃。

弱化这样一个噪音的这样一个方案，比如说我们有这个vehicle shape design，这个vehicle就包括了这个大道，飞机，小到一辆车，它的外形都可以被优化。

通过这样一个降低升压的这样一个方法来优化，来得到我们最好的这样一个降噪，还包括呢我们可以做影院的一个设计，在这个过程中呢，可以通过这样一个升压的分析和计算，去找出最好的。

在这个影院中构造这个声音的这样一个回响，这样一种方案，以至于我们可以做这个potential flow对吧，我们potential flow方面，我们就有有这个飞机的记忆的这样一个设计啊。

有我们的这个呃波浪的这样一个计算对吧，以及呢我们在湍流方面，我们有这个vertex method呃，那好吧，今天的这个分享就到这里了，哇真的非常棒，我觉得我今天听了一下。

我又对这个fm和p p m理解更深入了一些，之前看师兄好多西瓜和论文，感觉里面那些给每个每篇效果都非常的棒啊，包括其实从其实我记得你博士期间是三篇论文，其实都是哦。

有可能有两篇是vertex method对吧，我记得第一篇，第二篇，第三篇是一个呃是一个big flip的一个一个paper，对对对对对哇，我觉得还是每一篇效果都非常的惊艳。

我记得我第一次看到p p p m，大概是我大二还是大三的时候，那个时候感觉对各种东西也不太理解，然后这一转眼都可能有哇，这一拳也好多年了，好多年了，有同学问哦，窝方法如何处理固体边界哦。

五个方法处理固体边界的时候是这样的，求解边界积分方程，就是说我看我这儿找不找得到呃。

我这没有讲，但是我们回到这个边界园的那一个slides哈。

我回到边界源的这个slides，其实在窝方法处理边界的时候，我们是这样的，我们希望找到一个在这个边界上的，这样一个窝分布，或者这样一个电子分布，这个区别于这是一个势流，你是要做这个呃。

你需要做粘性流还是非粘性流，我先说这个无粘流，无粘流的话，其实我们只要保证法线方向是不穿透的，就可以了对吧，法线方向是不穿透的就行了，所以这个时候，其实我们需要找到一个什么条件呢。

我们是想在这里找出一个电子，你可以想象它就像在这里找出一个电子的分布，使得这些电子呢它可以introduce一个速度场，一个试流场，一个potential flow的这样一个厂。

这个potential flow的这个厂的刚好可以去干嘛，刚好可以去cancel掉，刚好可以去cancel掉，这个边界处的这些，本来会这个穿透的这些速度对吧，所以你立的方程呢。

就是这些东西产生了一个边界速流场，然后你让他去乘以这个法线方向ni对吧，它正好可以等于负的对吧，负的现在你已知的这个呃速度点乘n对吧，这个就是你要呃，你你对于这个无电流要解的一个方程，那对于粘性流呢。

你解一个额外的方程，就是说你希望在这里再找到一个这个沿着这条，沿着这个边界的butt sheet，沿着这个边界的一个vtec，使得它刚好可以in roduce一个速度场，能够去干嘛呢。

能够去cancel tangent方向上的分量，能够去cancel tangent的方向，那这个就是你找出来的这个这个边界，出的这个值了，那么这个边界处的值呢，接下来你再通过这个热扩散的那样一个方程。

或者说呃那个叫做什么，在这里，这因为在一个纯粒子的方法中，你是不太好不太好做热cosm，所以你可能要通过一些randow，或者通过一些呃这个能量总和不变的啊，这样一些扩散的方式。

把这个上面求出来的这个vertex系统，在这样一个强度，给emit到周围的这个涡流里面去，那这样的话呢你就能够形成卡门涡街了，对ok，所以相当于是你还是呃在你得在固体，比如说你说卡门国界的话。

你得在这个圆柱的内部去给他先求出一个呃，你就假设它没有固体，然后你求出一个就先假设没有布局什么分布，然后你就你这样吧，你就是这对我你你n时刻，假设你现在在n时刻，你所有的涡粒子已知对吧。

所有的速度已知对吧，这时候你先做什么呢，你先做对流对吧，你把握粒子走一下，然后你再去干嘛，你再去从涡粒子里面求出一个速度对吧，那这个速度呢，它在边界上不一定是严格满足的对吧。

或者说一定不是严格满足的对吧，那么它就有一个分量要被你cancel掉对吧，要看收掉多少呢，就是你去把它点成tangent或者点成法线对吧，嗯就是你要看售掉的量，那这个量由谁来cancel呢对吧。

由谁来cancel，那它就由谁，就由这个边界处的电子引发的一个视流场，或者说在边界上的这样一个vtex sheet，引发的这样一个速度差来cancel对吧，i c当你求出这两个东西了以后。

那个电子厂你不用管对吧，他就是个四流场，它无限的剩下的那个vertex 呢，你再通过一些热扩散的形式，把它上面的这些那个，把它上面的这些窝料给扩散到你的系统里去。

你就你就你你就已meet了新的涡粒子了吗，这边有一个边界条件，它是vertex的source，对对对，他就是一个探索，在他的边界上，对啊，就是一个butt s，假设这是一个无法移边界条件。

他会搓出这些word works，对对ok，有意思有意思，对这个这个我觉得可以这个写起来应该还好，如果其实就是一个找一个普通的vertex solver，然后在边界上不断去see的这个。

以某种方式去c的这个窝，这个窝的例子应该也可以搞出一个卡门，我姐在对，事实上对于二维的这个圆球不对，二维的spr和一个给定的那个速度啊，这个东西是有个解析的表达式的，你得找一下可以找到。

然后你就按照这个去去seed就行了，然后出来就是个卡门涡街，对哈哈，ok而且那个卡门过街特别的炫，就是各种湍流，各种turbent，特别的特别的是就对对对，是之前之前。

我们作业林里面有一个同学写了一个卡门，我借用lbm写的，然后还被中科院物理所啦去写了一篇科普文章，还是很有意思嗯，ok那么如果大家暂时还没有问题的话，大家可以在微信群里面和信息师傅保持联系啊。

我感觉这个sm还有这个pv p m啊。

这些快速学习方法还是挺有难度的，我觉得他的idea应该还好。

和mod其实差不多，但是你如果真的要实现它还是很难，理想情况下，我应该写一个这个，但是我看到他有什么m t m m t l l t l l t。

p p two p这些加起来实在是太麻烦了，而且很多各种常数啊，这些百米可以分享100，太会写，对三维的就更恶心，三维的你要去看一下那个公司都想上吊，那个那些在呃。

那spherical harmonic expansion去展开了，我靠那还要做rotation，就是三维在做l to l，还是做m two l那个transform的时候，为了加快它的一个计算量啊。

两个spherical harmonic expansion，先做一个对齐，先做一个那个旋转，使它对齐到半径方向，使它对齐到半径，然后再去做一个那个转换，就会那个系数就会更快，不然的话会慢一点。

就有很多像这样的那个很恶心，我觉得所以我说我是真的觉得差值应该是个好。

好，对对对，就是插值方法说不定会让他减少很多，但是m t m这些还是得做的呀。

因为很多是m two l就不做了，m two l不做了，只要做m two m就行，意思就是m t m还得做l dl不用做了，对对对对对。

ok那那我们就再次感谢一下星星师兄。

给我们分享的这个内容，我觉得真的非常的质量非常高。

而且非常硬核，从这个手写的pp就可以看出来。

真的是在百忙之中给我们来抽空讲这么一讲，真的非常不容易。

再次感谢感谢，ok那么我们就能保持以后的沟通嗯。

GAMES201：高级物理引擎实战指南2020 - P6：Lecture 6 线性弹性有限元与拓扑优化 - GAMES-Webinar - BV1ZK411H7Hc

ok那么我们就开始好，还是欢迎大家回到我们的课程，我们把电脑升级了一下，应该不会再有卡顿的问题了，然后还是希望大家多来听听直播，因为及时提出问题嘛，然后我们收到反馈，这个我收到反馈以后也会调整难度。

这样大家讲的讲的大家也容易理解一些，那么今天我们来讲一讲，这个线性有限元和拓扑优化，其实主要是讲线性有限元，拓扑优化呢只是顺带提一下，然后课程的附件里面有一个代码包，代码包里面呢它有这个嗯。

线性有限元的呃，这个top优化的一个代码示例，然后大家可以跑一跑，还有一个手把手教你做特殊优化的一个pdf，那个是mit的一个计算机图形学的一个作业，然后当时我帮他们写的这个代码，然后借这次就搬过来了。

给大家可以手把手的去做一做，当然这些作业我们也不批改批改，大家都不用提交，可以只是你有时间的话，你可以做一做。

没有时间的话，你可以看一看嗯，意思一下就可以，ok那么作为一的获奖选手呢，我们已经在植物上面公布了，大家可以去看一看优秀呃，作品的优秀同学的这个作业，然后和以前一样，代码都是直接。

一般get close一下就可以装菜就可以运行的，所以大家可以看一看，跑一跑呃，看看别人的作业是优秀的同学是怎么做的，然后这样自己也可以学习学习，那么当然了，肯定有很多同学说这个时间太紧了。

我没有赶上deadline，或者最近有事没时间做，这个完全不用担心哦，因为我们呃你加把劲，一个月以后就可以把它作为作业二来提交，作业，二呢又回到这个开放作业的模式，和作业领差不多呃。

只要是一个物理模拟器就可以，只要我们没有标准标准，就是你自己满意就可以，然后有很多很多个选项，你可以实现一个可交互的2d的物理模拟游戏，你你也可以去专门优化性能呃，你可以提高粒子数啊，或者提高网格精度。

因为高性能的物理模拟器，当你粒子数多了以后，网格比较细的以后是很多时候呃，很多是很多低性能的物理模拟器的东西，就不能再用了，高性能物理模拟器啊，里面有很多东西呃，并不能把地形的照搬。

比如说你本来是用稠密数组存在你的线性系统，但是你网格多了以后，你就可能重名速度变得非常大，内存存不下了对吧，可能啊你就希望用一个稀疏数组，这个时候呢就会意识到，你要做一个高精度的物理模拟，你要做挑战。

你要面临的挑战其实是非常非常多的，除了性能，你也可以去实现一些高精度的格式，比如说你以前是用的sama rg，你现在可以把它换成比如说bf c c，或者买comic这个高精度的格式。

或者你也可以用一些election，reflection scheme之类的去降低数值点性，使得你的烟雾模拟看起来更加的呃怎么说呃，充满活力一些就叫energetic，这个这个让我想起来一个笑话。

因为呃很多时候计算机图形学里面，大家有的时候你做这个sober，你的目的是去迎合你的导演的需求嘛，所以呃我想起来就是kaths是他是dream，以前在dream works后来到了维塔。

然后他们是负责vfx，他在那边就他有一次给大家介绍它的数据结构，叫vd b，然后他就讲到一个笑话说呃在统计学里面，这个很多时候大家希望艺术家可控性就是呃，instead of。

你去实现一个特别高精度的火焰，有的时候导演会让你实现什么，导演会让你实现一个生气的火焰好，这个时候你就要想想呃，怎么样去提高你的这个呃simulator，怎么样艺术家去控制你simulator。

因为simulator的它的背景，它的底层是partial differential equations，实是这些偏微方程组对吧，但是有的时候艺术家并不会直接去操纵偏微方。

他可能会想用一些其他的办法去操纵这个火焰，所以这个艺术可控性，还是一个很重要的一个事情，但是这个也可以作为选项四，然后选项五呢你可以自己想想，有什么别的有意思的方向，可以自己去研究研究。

或者找一篇比较好的syrah的模拟方面论文，然后自己去实现一下也可以，当然你可以基于自己的或者别人的作业一，然后呢，我们的评分标准呢是以新增的这个原创部分，然后还是建议大家组队，因为组队的话呢。

首先你如果组队的话，你如果拿了奖，那我们会你队里面有多少人，我就会给你寄多少个杯子，所以呃记多少个课程纪念品，所以这个组队的话，从得到纪念品的效率来说是比较高的，但是这个是一个呃不是最主要的。

最主要的是你组队了以后，你会在做这个写代码的过程中啊，你会呃提高和别人合作这个能力，这个还是挺重要的，然后这个截止日期呢是北京时间8月15日，23：59，也就是说接近一个月以后，那么接下来呢。

我们基本上就进入了课程的后半部分，前半部分呢大家嗯应该来说还是有一些难度的，很多同学反映，特别是推这个影视时间积分器的时候，这个头都秃了，难度还是挺大的，这个其实非常正常，大家不要受这个打击。

我一开始推一个时间分期的时候，我一我一开始推，我记得是大三的时候，我去实现，大二的时候去实现一个m p m的影视，实现积分器，那个真的推的我是头秃，但他经常是什么呃，四阶张量之之类的。

推起来特别容易推错，然后推了好久才把他搞对，然后弹簧质点其实也挺复杂的，和m比起来他更不规则，它会有一些对于什么，比如说呃x取了normalization以后的向量，关于x求导之类的这些项目。

这些项一不小心就会推错，然后这就是为什么之前也推荐大家组队，因为组队的时候，这个公式大家可以互相检查对吧，这样实际上两个人都退缩的概率，就会比一个人退缩概率要低很多，那么接下来今天我们讲有田园。

然后接下来两讲呢，我们会讲混合欧拉拉格朗日视角，这个其实包括一些常见的方法，比如说particle excel或者flip或者iphone，particle insl，然后也会讲到物质点法。

物质点法应该是这两讲里面的重头戏，到了第九讲，我们会讲一讲高新计算与物领情，其实就是一些常见的优化物理引擎性能的技巧，那么8月10号呢我们还是像之前一样，我们空一次，大家留一点时间。

可以消化一下之前讲课的内容，也给大家留一些时间完善自己的物理引擎，最后一讲呢我们进行总结，可能会点评一下优秀作业，可能会回顾一下整个课程内容，最后可能会讨论一下呃，物理物理里面一些还没有解决的问题。

ok那么这个是课程安排部分，接下来呢我们就到了今天的重头戏啊，有些部分我还是把这个图标拖到左下角，这样不会挡住，好开奖之前先喝口水，这个后面肯定会讲的口干舌燥，因为今天其实公式稍微有点多。

那么这部分呢我们分为四个内容，首先我们大概看一看有限元是干啥的，然后后面第二第三部分呢，我们举两个简单的用有限元来离散化的例子呃，首先我们从大家最熟悉的泊松方程，标量泊松方程来开始讲，然后呢。

我们再介绍一个柯西动量方程的优先原理散化，当然我们呃讲的是最简单的，就是说呃基本上只有一项的这种动量方程，没有为没有外力，没有这个呃质量的flash就是没有这个密度的变化，所以其实只有一个t。

然后最后呢我们大概介绍一下拓扑优化，那么有限元方法，它其实是加六经方法里面的一种，然后什么是加六经方法呢，加六经方法就是一种呃，把强行的pd转化为弱性的一种方法，嗯这个后面什么是强行，什么是弱性。

我们会再细讲，在有些园里面的，连续的p d e会一步一步经过变换，变为一套离散的线性系统，然后大家呃就可以用最喜欢的线性系统修剪器，像比如说工作梯度呀，mtv的呀，去sop这个有限元的系统。

那么经过计算机图形学里面的简化，一般这个有限元他还是就这么几个套路，可以看到这边有1~5步对吧，首先我们把强行的pd转化为弱转化为弱性，用一个test function w。

后面我们会提到什么是强行的pd呢，强行的pd就是每一点都成立的pd，它是一个连续的形式，它在空间里面连续的每一个点都成立，那么这个大家可以感觉到他肯定是呃，不太在计算机里面是不太容易实现的对吧。

因为你计算机里面只有有限的呃，存储和有限的计算，所以你是没有办法去做这么个连续的表示，呃然后刚才说的第一步，然后呃哦对我还没说什么是弱形若行，就是说我们把这个强行pd它的残差项。

用一个test function加权，然后积分这个后面我们会提到，那么第二步呢我们做这个分部积分，分部积分以后，我们可以做一些实际上把系统进行一些简化，我们一般可以把呃get it的去去到里面。

pd里面的二阶导数项，然后第三步我们用一个散度定理，然后再进一步简化纤维方程，同时呢在这一步里面，我们可以去得到第二类边界条件，就是诺伊曼边界条件的一个显示的形式，然后呢我们在做离散化，离子化。

也就是说我们把这个空间里面连续的场，用一系列离散的点来表示，那么大家应该已经很熟悉了，线性差值对吧，双线性插值，三性差值，其实这些都是非常常见的呃，用离散的量表示连续量的方法。

有些园里面有很多很多种插值的方法，但有些园里面一般叫basic function，它有很多种basic function可以选择在这边这一讲里面，我们其实只考虑tt这个线性的basis方式。

后面会有一些图，大家可以有一些直观的理解，那么你把这个连续的常用这些呃，离散的basis function的加权组合来表示，以后你就得到一个离散的系统，那你最后可以得到一个大家一般会常常见的。

叫做stephens matrix，但是他这个不一定是stephens，在很多时候它可能是一些其他的呃，物理意义，并不是总是说是一个刚性这么一个东西。

然后除了这个matrix才会得到一个right hand side，就是这个load bor，最后一步我们求解信息系统dj f m非常重要，因为很多其他离散化方法，像有这个叫物质点法呀。

他其实也是一种加勒定方法，他的思想和m和fm其实非常非常像，只不过在npm里面呢，它没有显示的元素，它是一种叫element free grian，它是叫呃无元素加到金发，ok。

那么我们还是从大家最熟悉的二维泊松方程，来看解，那么普通方程是啥，大家应该还记得吧，就是呃我把这个我有一个标量场，这个u我给他求一次梯度，再求一次散步，然后等于零在整个空间里面。

那么这条空间呢它其实还有四条边对吧，那他在这个问题里面，我们假设前它的这个三条边，它是狄利克雷边界，什么是狄利克雷边界呢，底力科技可类边界，就是说我直接指定这个u在这边界上，这一点的值。

我用一个fx去表示它的边界上的值，那么下面一条边呢，它是一个诺伊曼边界条件，也就是说呃我们不直接指定ux值，而是我们指定ux它的导数的值，当然我这边得和边界上的法向量做一次点。

成了我们gx表示它的洛伊曼边界条件的值，一般大家也会把德雷克雷边界条件，叫做第一类边界条件，然后诺伊曼边界条件呢叫做第二类边界条件，然后我们这还是一个连续的问题，左边这个是连续的问题啊。

然后我们一般来说呢，在统计学里面，很有可能会用一些网格去把它分成很多，很多个这个小的部分，然后同学里面大家比较喜欢用regular的，这个叫规则的四面体四边形网格，四边形网格，然后很多同学会问。

如果说我把这个100和u11 连一条边，那我会得到两个三角形，就得到一个三角形网格，那这个可不可以呢，这个也完全可以，实际上很多时候大家会更喜欢用三角形网格，三角形网格有很多的优势，比如说它可以呃。

你可以很容易地做网格的密度的渐变对吧，就adaptivity，然后三角形网格也可以，这边我其实展示的是一个方形的边界，如果边界是一个圆形或者是一个其他形状的话，你可能会更想用三角形网格。

因为三角形网格能够更好的贴合你的边界，当然我们今天这一讲，我们所有的数涉及到数字的部分，我们都是用正方形的元素来表示的，ok那你可以看到我这边其实分了多少个元素，分了3x3=9元素对吧。

那多少个自由度呢，其实我们的自由度是存在，每一个元素的四个顶点上，那我们其实有4x4=16个自由度，也就是说16个浮点数来表示这个呃，这个场具体怎么表示，我们后面会提到这有啥用，这其实有很多用。

光是一个标量的泊松方程，你就你就可以用来解流体模拟里面的压力投影，对吧，这是一个最常见的用处，当然一般大家在流体模里面用fm用的比较少，一般用图形学里面常用的是final difference。

是有限差分和final morning是这个有限体积，好刚才我们介绍这个问题，那么接下来我们来看看到底怎么离散化它，首先我们得把这个强形式给它，转换成一个弱形式，怎么个转换法呢，我们想象呃。

我们有一个任意的二维标量函数叫w x，这个w x呃，我们把它和之前的强形式去乘起来，然后做一个积分，然后在我们的整个东北ea上面做一个积分，那么我们有这么个关系式，这个关系式就是说。

如果说我们这个pd在每一点都成立，当且仅当对于任意的test function w，它这个积分都等于零，那么我们来直观上理解一下，如果我们要从左边推到右边，其实很好理解对吧。

嗯这个如果说左边这个项永远都是零，那你这边是w不管这个w不管是什么，你乘上一个零都是零，然后你再和da去再积分得到还是零，所以从左到右是非常tribute的，但是从右到左的，那么这个就需要一些想象力了。

如果我们存在一个点，它的呃强行是这个三，他这个plan operator不是零，那么怎么办，我们可以总是可以，想方设法构造出来一个test function w，使得右边的这个不等于零对吧。

然后我们这个叫做利弗利，把它换一换，就可以从右边推到左边，当然这个构造这个只是一个直观上的构造，你真的要去从数学分析上去证明，它还是要下一番功夫的，ok那么我们接下来要做一个什么事情。

我们接下来希望把这个强行式里面的呃，二阶导数给它去掉，用我们的这个test function咋去呢，我们来看一看呃，首先我们要用一个技巧叫做分部积分，这个大家应该都非常熟悉了。

其实或者你不调到分部积分也可以，你可以叫他乘积的导数对吧，你看这个右边它就是w乘上grab u，你要算它的散度，那你这边呃求导数的话，其实就是两个term对吧，你先给w求导数，然后再点乘上u的导数。

或者呢w乘上这个这个u的这个二阶导数，或者说他的npsl trader，那么你看我们做了这个变换以后有什么好处，我们得到了一个这一项，这一项，他你看上面这个方程这一项它其实是零啊，那你就可以做一步简化。

你把这一项给他踢了，那你就得到这么个整式，那么你会发现吗，上面这个公式等于零，当且仅当这个公式二是成立的，那么接下来我们把test function，我们不apply在原来的这个上面，这个公式上面。

我们把它apply在公式二上面，就可以得到公式三对吧，公式三你可以看到呃，uda plant的一零，当且仅当呢，对于任意的w我们有这么个公式成立，有人会问你这不还是有二阶导数吗。

你这个地方不还是有一个二阶导数，后面我们会讲这个地方可以怎么把它去掉，ok那么公式三到这边来变成公式四，我们怎么把右边的这个二级老师给他去了呢，其实也很简单对吧，你看到它其实是一个散步的形式。

那么一个体积里面散度的机会等于什么，等于它边界上的啊，这个向量函数在呃和法向量点乘，绕这个边积分一圈对吧，这边其实用了一个散度定理，那冥龙彩龙定理以后呢，你会发现我如果在二维里面。

那么我左边的其实就是一个呃面积分，右边其实就是一个沿着它的边界，一个线积分好，那么我们所有刚才讨论的，全部都是局限在一个呃连续的空间里面，接下来呢我们要做一个事情，就是说我们用一系列的离散的函数。

去表示这个连续的函数，那么这些离散函数呃，我们怎么来用一些离散的basis方，也就是说这个基函数来表示一个离散的呃，来表示一个这个原来这个u呢其实也非常容易。

那么首先我们有一系列的呃base function，假设我们其实是有四个道，假设我们是five 0123，那么假设我们的定义域啊，就是从这个原点到这一点，那你其实可以看到我们的f0 和f3 。

它是只有右半部分或者左半部分，那个fi一和i2 ，它是有两边都有，它，这个我们这边用的是一个，最最最简单的线性的计函数，你可以看到公式六，其实直观上的理解也很容易了。

但其实就是把我们的这些奇函数根据优先加权，或者说scale一下，然后我们再求个和，那求个和以后，你会发现我们其实就得到了一个呃连续的呃，也不叫连续的，这个叫peace by senior。

是分段线性的这么一个u的函数，它其实在很多时候能够做的挺好的，挺好的，近似的，那么在这门课里面，我们专注在这个线性的计算中，上面，有同学说我这边能不能用一些更加fancy的一些，我不我不喜欢直线。

我把它换成曲线也行不行，然后完全可以，那么那种情况下呢，那你其实就整个系统就会变得更加复杂一些，那有些同学会问了，那你说这么复杂，这东西不就是线性差值吗，对他确实就是线性差值。

你如果嗯认为这个幽灵是就是零点的这个值，因为就是一这一点的值，那它之间的值啊，确实就是现性差值出来的，那一维里面是线性差值，大家可以推导出来，二维里面是什么，就是双线性插值，三维里面就是三线性插值。

在统计学里面，大家其实一般情况下呃linear或者by linear，trilinear，差值一般就够了，那么有的时候会有些呃文章里面会用高阶的，那也可以，但是相对来说比较少，那刚才看的是一维的情况。

到了二维情况会怎么办呢，注意这个是一个二维的图啊，差不是三维的，我这边啊这个图里面是用高度来表示函数的值，但是这个图是偷的，这个图是从cal的网站上面找过来的图，我实在是没有精力去画这个图。

这个画图实在太费时间了，大家可以呃，如果这门课结束以后，大家可以去看一看council的这篇介绍fm的文章，他这个其实还是写的非常好的，当然人家是商业fm软件，他写这个文章一方面是介绍f1 m。

一般一方面是也是宣传自己的软件呃，大家可以看看他这个里面图都非常精美，ok那么二维的basic function长啥样，如果你是三角形网格的话，那你basic function大概就是长这样。

它其实是一个大家一般叫它一个帐篷函数，但非常像一个帐篷，它每一面都是三角形对吧，然后它在中间高，然后相邻的道上面，相邻的这个顶点上面全部都是零，那么呃你会发现啊，他其实这个杯子方式。

他其实大部分basic function是不相交的，刚才我们看的这个例子，它是相交的，把他这个大i m fig它是呃有一分是重合的，但是大部分情况下但是不重合，这个有什么意义。

这个后面你就会发现它不重合的话，就对应的一个是你的线线性系统里面，这一个啊零这么一个值，那么他如果说重合的话就是非零，那么既然大部分basis function它是不重合的。

那它说明这个线性系统里面基本上大部分都是，所以说你得到了建议系统啊，它其实是一个稀疏的线性系统，ok，那么我们还是回到我们最喜欢的这个，四边形的元素，这边我们其实用的是工工整整的。

方方正正的这个长方体的元素，这边就要用一下你的想象来去visualize这个basic function，这边其实在二维里面，它其实就是一个双线性插值的base function。

或者说它就是一个呃一个曲面版的tt function，tfunction，它是一个帐篷，但是它的它会有四个面，对它的每一个面是一个曲面，为什么是曲面，而不是而不是平面吗，刚才我们看的这个三角形元素。

它不是平面吗，ok那么，这个其实是涉及到blinear或者trainer，interprelation里面有一点让人迷惑的一点，就是blinear interpolation，虽然说它叫blinear。

但是它是关于x是0。2，关于我是零零，但是它关于x y它就是二次函数，所以你得到的并不是一个完全0。2的对吧，它其实是因为它是把x方向的这个linear，和y方向的linear给它乘起来。

所以得到了本质上还是一个二次函数，然后总体来说啊，那一般你肯定是呃，比如说u11 这一点，它的形函数，那它肯定是和刚才一样，它在周围的这八个点上面都是零，中间一个点是一。

啊对我刚才可能用了一下这个新函数，这个菜其实就是或者或者叫基函数，我一般可能比较喜欢用g函数，因为basis function它一到中文应该叫基函数，它其实呃但是有些地方它会有混用不同的称。

有些园里面它这个地方啊，有些地方它的这个术语它不是很规范，有些人会有一些强迫症，他必须这个地方用shift方式，那个地方用business function啊，但是我们就不去做明显的区分。

ok那么接下来我们刚才说了这个ux啊，我们其实可以把它表示成若干个，basic方面或者基函数的加权求和对吧，那么我们再回忆一下，我们之前想soft这个方程四，其实到这边变成方程八了。

那么其实我们想求解一个u使得什么呢，使得我们这个方式能够满足对吧，右边怎么来的，还记得吗，右边是散入定理搞来的，左边是怎么来，左边是分部积分来的，ok那么我们简单的substitution。

我们简单的把这个ux带到，把这个离散码ux替换到呃方程八里面去，就是说把七塞到八里面就会得到一个九对吧，这个九你会发现它其实就是说呃，当然我们只换左半边，右半边，这个因为我们是不换的。

因为右半边我们后面有别的处理方法，ok那么，这个里面其实还有一个w，我们一般来说w你也可以用这个base function来表示，但是还记得吗，我们这个w是任意的w，然后我们这边只考虑它的离散的情况。

就是说这个w也是由basic方式来表示的，那么你这边如果要取任意aw，其实啊这边相当于说你只要取任意的basic function，作为w就可以了。

这边其实说呃因为w本身也是basic function，给它线性组合出来的嘛，所以我们只要取它的g就可以了，ok呃这个地方可能得稍微想一下，但是呃我们可以非常安全的，把任意w换成任意i呃。

人鱼发来这个地方得想一下，但是还是比较容易想通的，因为你如果不取极限g函数的话，反正它也是线性相关的，你只要取线一向线性这个无关组就可以了，那么我们把w换成fire了以后，我们就可以得到这个公式十。

好公式十到这儿变成了公式11，那么我们把呃接下来做一部非常有意思的事情，我们需要把这个ug啊给它，从积分里面给它提取出来，提出来以后，这其实就是一个比较脆弱的交换顺序对吧，那么我们得到的公式12。

我们只换左边这部分，右边我们还是保持它原样，那么你会看到公式12，变成一个我们非常熟悉的形式，这边未知数其实是u g，就是呃这一点你要求的压强的值，或者一些什么其他东西的值，然后。

呃你看ug左边这么一整坨这么一个括号，它其实啊和尤其是没有关系的，他其实就是一个如果知道的话，知道了，这其实就是一个常数项对吧，那么右边呢右边这个我们会后面我们会讲，它其实是一个啊。

可以用诺伊曼边界条件给他推出来，ok你如果把它写成矩阵形式，你就会变成这个ku等于f，那ku a等于f是其实就是ax等于b，在有限元里面的一个大家喜欢用了一个别名，它其实就是ax比它就是一个线性系统。

那么在k有限园里面，这个ku f都有自己的名字，k一般大家叫做这个stiffness matrix，或者叫刚度矩阵，然后u呢u一般大家叫做嗯自由度，或者叫做解向量，f呢叫做load bor。

就是说嗯这个其实更像是在固体，在解这个固体小形变里面会用到的这些术语，但是这边我们呃大家习惯了，所以我们就照用，还记住这几个名字，ku f，那么还记得我们刚才前一张size里面有这个。

有这么括号里面这么一坨，那么我们它其实就是下面这个矩阵形式的，k i j对吧，那么k i j它是啥呢，它是呃在整个冬梅里面积分，他这个fi和fg他们的梯度点乘的积分。

由于我们这边用一个简单的base function，所以其实他这个积分不难计算，那么在一些更难的basic function里面，大家会用高in queror，但是这边我们简化以后。

我们其实手动积分一下也没完全没有关系，那么假设我们中间这一点是i，那么fi其实是什么，大家可以想象一下是一个帐篷函数，对吧，那么他的这个帐篷函数啊，如果说我的j和i是一样的话。

那其实就是两个相同的账目函数，在这里面算这个fia呃和fj，然后算他们的梯度，然后变成然后积分对吧，那其实由于我们的basic function它的支撑是有限的，它只是只在一个很小的范围里面，有。

他这个base function只在周围四格里面是非零的，在更远的地方全部都是零，所以我们其实只要考虑什么，只要考虑他周围的四个元素的面积，里面的积分就完全ok了，ok然后如果是i等于j。

你会得到一个这个是8/3，那你如果j在这边还在这的话，你会算出来的是负呃，1/3对吧，然后周围全算出来全部都是-1/3，当然这里面还有一个一些像比如说dx啊，这边我们就忽略掉了。

你如果说把它再scale一下，就可以得到一个九点的，拉普拉斯的aplan operator的一个dl，中间是八，周围全是-1，那么还记得我们之前曾经用有限积分，这个有限差分。

final difference来呃，推出过一个五点的拉普拉斯stl吗，就是-4，然后四个角是一对吧，这是呃这可能有正负号的一些问题，然后但是大家一般也不care。

然后你会发现它其实同一个operator，你用不同的方式可以算化的，可能得到不同的cs，当然其实大家呃在我之前提到图形学里面，一般大家做这个流体，还是更多时候是用有限差分或者有线体。

所以用这个stl用的多一些，前面这个中间是八，周围是一的那个三层，反而用的比较少，但是呃播送方程作为推有限元的例子，还是非常非常好，因为它真的是最简单的，不能再简单了，还记得我们之前一直有一个。

像我们一直没有去讨论的，就是这个方程的右转向额，我们用散步令，你把这个散步的积分换成了边界上的，这个应该叫做啊，我忘了这个叫什么名字，反正是边界上的这么一个梯度的积分对吧。

还记得我们这个呃u的梯度点乘上normal，这个是一个什么项，这其实就是我们之前提到的这个，诺伊曼边界条件对吧，所以我们其实在这边，在这个地方我们就解决一下边界条件的问题。

首先我们考虑一下这个dva类边界条件，也就是说叫第一类边界条件，那么呃其实非常容易，这个其实enforce第一类边界条件非常容易对吧，我们只要把相关的ui，也就是说在边界上的那些自由度，它的值设置成嗯。

这个第一类边界条件给定了，f在xi xi就是ui那个点的坐标，我们直接设置，为什么呢，因为你如果是，你如果用线性插值的方式来理解的话，那其实确实那一点上的值就应该是f的值。

那么诺伊曼边界条件稍微稍微复杂一些，你会看到我们之前是用gx来给定的呃，user梯度点乘上normal对吧，那么你会发现这边其实这个项目，你也可以直接塞进去，那你把这个g塞到这个14的右端项。

你可以得到一个呃，最后你会得到一个呢，你low vor里面就会有一个非零的一个entry，所以其实你会发现，其实诺伊曼边界可能会导致你丢单项，有这个非零，在于有些有些园里面，大家会在这个。

固体里面大家又来向一般会叫做这个traction，好像应该叫做牵引力，就是说实际上就是你这个塞给他们力，所以有的时候大家会说，我推一个弹性的一个材料，我其实给他apply了一个诺伊曼别的条件。

这可能会有人让让人摸不到头脑，他明明就是一个呃推力，怎么就变成了诺伊曼边界条件呢，啊他其实从这个地方就可以看出来，好那么接下来我们进入第三部分，刚才我们讲的这个是一个泊松的方程。

我看一看有没有图像有问题，啊对，刚才说这个traction有有同学说叫做逆变焦，或者或者叫做面离啊，这个嗯还是非常好，这个有准确的术语嗯，lovector对。

lovector中文好像叫荷载荷荷荷载荷荷载啊，两个都是多音字，我就乱读了，应该叫荷载啊，反正呃就这样了，嗯哦对再合啊，anyway，有能不能把这个拼音的声调也帮我打打出来，这两个字都是多音字。

ok嗯对，有同学说在算力的时候，在力学里面第一类边界条件是什么呢，第一类边界条件就是这个点的位移，如果说你会经常说我这个固体这个地方，我去定了个钉子，把它嗯固定在了平面上，这就是第一类平台的点。

然后我如果呃给加一个弹簧在那去拖拖沓的，这个就是第二类边界条件，对这个我们后面其实会提到，好，那么我们接下来讲一个稍微稍微复杂一点的，有限元离散发，这个就是呃线性弹性的这个有限元。

这个其实呃是从什么推出来的呢，是从这个柯西动量方程推出来的，那么当然上面这个方程稍微有点复杂，你可以看到其实有一个它也很容易理解的，还是左边的是一个v的拉格朗的导数对吧，就是说你的加速度等于什么呢。

其实就是f等于ma对吧，呃或者或者说这个其实是呃a等于f除以m，那么右边像是什么，右边是嗯就一个16/6分之一，其实就是f等于ma里面这个mg一项，然后这个地方有一个稍微有点难理解的。

就是说我这边要算呃柯西stress tensor的散度，这个可以sensor，它是一个3x3的一个张量，或者说叫或者是矩阵了，然后它是对称的，在三维里面就是3x3，二维里面就是2x2。

然后算tsa的二阶tensor的这个散度，后面我们会提到有没有什么容易理解的方法，然后再加上一个g g是什么，g就是体积力或者说重力加速度，ok那么我们先把他大刀阔斧的简化一把。

首先我们考虑的是这个准静态的过程，就是说它的速度它是零对吧，然后或者说就是一个微小形变，就是根本就没有和移动，它只有一个微小的形变，所以他v 10，说到这个准镜台，我想起来一个笑话。

我有一次我本科的时候学计算机，然后呃我当时在摇班，然后我隔壁的同学，因为姚班有学计算机或者学数学，或者学物理的同学，有一次我起床，可能是前一天写代码写的时间太长了，然后起床就落枕，脖子很疼。

然后呃我就跟我室友提到这个事情，我室友他是学物理的，然后他就然后我跟他说，我这个脖子啊，它这个疼痛感是随着我脖子转动的速度成正比，然后他上来就来了一句，那你以后可以准静态转动脖子，这个什么叫准静态。

就是说你虽然在转，但是你每时每刻速度，就是你就是说你很很慢的转你的脖子，ok啊这个就是一个关于cos static的一个笑话，我每次看到cos static，我都会想到本科时候这个这个经历。

然后呢刚才讲到速度是零对吧，然后我们的density我们认为它是不变的，那么我们可以把这1/6这个term给他去了，然后也没有gravity，也没有体积力，ok那么我们把这些全部去了以后。

就得到一个呃最简单的形式，就是说呃克星stress tensor的散度等于零，那么我们这边的自由度是什么呢，我们还是用u表示自由度，这边的自由度是这个displacement。

什么叫displacement，就是说我的呃这一个点在施加了外力的情况下，他得有一个形变对吧，它形变了以后，这个点它的形变的位置，和原来的位置就会有一个位移，这个位移就是displacement。

注意啊，之前我们讲嗯，之前我们讲这个最简单的三角形f e m时候，那个其实是一个比较山寨的f m嘛，然后他根本就不是用呃，我们那讲讲的时候，根本就没有提到他是怎么推出来的对吧。

呃那个时候我们是用了一个fifunction对吧，然后by方式呢它其实呃我们当时那个fifunction，它其实就是呃新的位置，关于旧的位置，关于material space的位置的一个映射。

但这边我们不用那个范了，我们这边是用的是新的位置减去旧的位置，得到了这么一个displacement，这边要注意一下，和之前用的不是一个东西，那么其实这边你会发现，这个方程里面好像没有出现。

我们的degree of freedom对吧，没有出现有，那么是为什么呢，因为这个sigma它其实是u的一个线性函数，那这个中间有一套非常复杂的映射，那么当然是一个线性映射，后面我们会提到有同学会问。

这个东西到底是一个login还是alin的，但其实他这个是一个呃infantasia的，infinitesimal的一个declaration，它是一个呃微小微小的形变，所以你他根本就没有移动。

所以你这边说他是loona或者lin的，这个其实没有太大太多意义，呃我这边还是提到这一点，因是因为很多搞graphic同学会觉得我用了网格，我就是欧拉的方法，但其实并不一定啊。

他这个其实并不是根据你用网格，还是用例子来区分它到底是leona还是order，ok那么刚才我们提到这个coaches，stress tensor divergence，这个其实可能不是很容易理解啊。

我们这边引入一套新的符号系统，叫做index notation，这个系统还是非常有用的，后面我们推pm的时候也会用到这套系统，那么什么是index htation呢，就是说我们用alpha beta。

gamma这一坨希腊字母来表示它的index，那这样做了以后，我们就可以统一的呃，把x y z3 个轴或者w有四个轴的，我们一般就考虑三维问题，所以其实这就只有呃x y z了。

那么我们就直接用x阿尔法来表示，这个阿尔法可以是零的时候呢，第零个轴其实就是x，如果二分等于一得到的就是y就是第一个轴，然后背上如果是二分等于二的，知道就是z我这边其实不应该吧。

呃x不应该用x贝塔x伽马，我这边应该是用x0 ，x1 x2 或者x1 x2 x3 ，取决于你是从零开始输出还是从一开始输出，呃，总之呢我们是用这个下标，统一的下标来表示空间的轴。

然后呢用一个逗号来表示口径导数，比如说啊刚才我们的这个coach stress stress test，如果说我们是他是呃，阿尔法贝塔这么一个component，关于gma这个轴求导的话。

我们就用这个sigma下标是阿尔法贝塔，然后逗号伽马，ok记得逗号后面的这个伽马表示求导，然后这么做了以后呢，我们就可以得到一个，把原来它是一个vectannotation。

我们可以把它变成一个index notation，这样看起来就比较清楚，因为说实话这个33tensor散度，一开始他是什么意思，其实你还是要想一下的，但是用index station它就很容易表示出来。

ok那么其实你可以看到，我们如果考虑左边，他这个这个向量方程的alpha component的话，你可以看到，那其实有几项我们可以明显的很容易就换，比如说这个v可以换成v2 阿尔法对吧。

这个g体积的这个body force，我们可以换成g阿尔法，然后他的内力，它的这个内部的力呢，这这一特就是呃koshi stress tensor，散步除以除以o，然后他其实你可以把它换成对于贝塔求和。

然后实际上是stress tensor，他的阿尔法贝塔component，component，关于贝塔求和呃，关于贝塔求导以后再求和，ok，这么一个东西，他就说如果你一个下标出现两次。

那就会有一个隐藏的求和好，但这边我们就不用呃隐藏求和，我们就直接把所有球和全部都用西格玛写出来，这样就更加清楚一些，因为出现了两次，这个不是很容易定义，比如说你这边阿尔法这边有两个阿尔法，有两个贝塔。

到底是两个头球和呢，还是只对其中一个求和对吧，这个就呃对于初学者来说不是很友好，我们这边就假设，我们把所有的这个求和全部写了出来，好那么接下来呢我们用有限元去比赛划一把，注意啊。

这个和刚才的泊松方式比起来，稍微稍微复杂一点点，为什么它复杂呢，因为刚才柏松他是他是一个标量的一个pd对吧，这边我们是一个十张的pd，而且呢刚才的柏松，我们这个u是直接出现在房车里面的。

但是在这个有限元里面的啊，在这个呃柯西动量方程里面，我们还有一个mapping，这个main是把你的u map到sigma对吧，然后接下来我们推导了现在sigma上面做。

然后后面我们再讲怎么把u给替换上去，all right，那么首先我们和刚才一样的套路，还记得第一步是什么吗，有些人他是加料金法对吧，首先我们给他上一个test function，给它沉下来呃。

当然这个test function要注意一下，它是一个标呃，这个矢量test function，它是从这个二维平面到一个二维向量，那么我们得去把它关于贝塔求和对吧，那么其实就得到一个这样一个方程。

那么为什么要上test function呢，上test function以后我们就有两个像，有两个像以后我们就可以分部积分了，或者用这个导数的乘法公式呃，这个乘法的导数公式对吧。

ok那么我们来看一看我们刚才得到的这一项，这一项呢它分到左边就是两项对吧，注意这边我们把它括起来，整个对贝塔对于x贝塔求的一次导出，那么求导数其实用index station看起来就更加清晰了。

你看我要不是给把这个贝塔塞给左边，要么给他塞给右边对吧，那么我们其实其实还是一个叫什么，这个口诀叫做d u b u d v加v d u对吧，我记得当时微积分的时候学过这个口诀，然后ok做了分部积分。

还是和刚才一样的套路对吧，这边这一项我们知道它是等于零了，所以我们可以把这个term给他删了，删了以后呢，就只有右边的这两个项目，那就得到这个这么个方程，什么方程，那么大家应该可以反映出来。

这个地方我们接下来要做什么，接下来就是给他做积分，做积分以后，这边可以用散落定理，散度定理以后，就彻底的把它换成了一个面上积分，或者线上积分，左边呢就是一个面上积分或者体上的积分，ok。

那么接下来我们还是做一样的事情，有了公式15，到这儿变成公式16，然后呃公式17是和刚才一样的套路对吧，一模一样，就是说我们把连续的方程，连续的这个u这个函数把它替换成离散的。

用basic function给它支撑起来的这么个函数，ok然后呢呃到了18里面，我们就还是和刚才一样的套路，我们可以把这个w什么test function，给它换成basis function对吧。

我们不需要去测所有的w，我们只要测比较有代表性的那么几个w也行，有代表性的w是啥呢，就是base function，然后u呢还是把它塞进来，塞进来的时候注意我们只换左半部分，不画右半部分。

右半部分这个地方它是刚才提到的，只是一个面临，ok那么其实左边呢我们其实呃你会发现左边，由于啊我们这个左边是一个关于sigma函数，它并不是一个关于sigma的方程，并不是一个关于u的方程对吧。

所以我们其实还不太好直接把它放进来，那么接下来我们要做一个事情，就是把这个sigma和u的关系给它确定下来，这样呢我们就可以把左边彻底的给他，换成一个关于u的方式，就可以把这个sigma给他删了。

那么如果把sigma去了以后，那么其实左边很容易又可以变成ko的vf对吧，刚才我们讲到，其实我们需要把这个sigma和u给他联系起来，对吧，这个u怎么一步步到这个sigma，那这边有一个string。

有一个从啊u到string tensor的过程，string tensor是什么，我们之前提到这个这门课里面，所有的string都可以把它换成deformation，它其实就是一个。

但是局部形变这么一个东西，当然displacement和呃这个这个，displacement也是measure形变的一个东西，在这边我们其实要给他再求一次导数，才能得到局部的这个形变率，获得这么个东西。

然后这个string tensor这边我们其实就呃很容易的，我们可以做一个定义，它其实就是呃u的梯度，注意这边又有这个向量矩阵的问题了，他这个是虽然是一个呃三维的向量，但是你求的梯度里。

它就变成一个变成一个什么，变成一个二维的，就变成一个33的矩阵或者一个二阶tensor，那么就得到一个这个就得到一个这个string test，那么string cher呢。

我们可以把它变成再经过一些变换，注意这个lambda和mu是什么，lambda和muse，是我们之前提到叫lame parameter对吧，然后呃在统计学里面还是非常非常常见的。

特别是你在呃模拟弹性的物体的时候，ok那你得到这个呃，这个是一个柯西的这个stress tensor，它等于什么，它等于lambda乘上sin test的trace乘上i。

然后再加上两倍的mu乘上这个呃street cancer，当然这个可能看起来非常的confusing了，我们可以用index dotation来进一步简化一下。

那其实就得到这个e的阿法贝塔component，等于什么，等于12/2分之一个u的alpha component，对贝塔求导，加上u的贝塔对阿尔法求导对吧，然后呃，后面我们也可以用类似的方法去表示一下。

然后这边我们有一个这个dota方式，这个调查方式其实就是阿法和贝塔相等的时候，等于一，不相等的时候等于零，然后呃用这个调查方式呢，其实对应的就是identity matrix，说了这么多啊。

其实这些这些都不重要，重要的是什么，重要的是这个sigma是ud linear function，那说了这么多啊，终于我们得去build这个linear system的举证了，刚才说到呀。

这个西格玛是u的linear function，然后我们又有公式19，ok，那么我们这个地方其实还是可以推出一个k u，等于f了，但是这个地方要稍微注意一点，这个我们公式的这个左边算k这个term。

你这个u要提出来要下一番功夫了，他不是说直接可以提出来，他中间还要经过stress string matrix和string displacement matrix，后面我们会提到这两步处理。

才能得到这个真正的k，所以真正的这个fm它的呃弹性fm里面，他的这个k啊，他这个矩阵是要把它组装起来，还是不是特别容易的，你要去考虑他这个senal size，就是说如果说我们没有边界条件的话，呃。

对于每一个u，它会和多少个其他的e u的component的接触呢，也就是说你这个k里面那一行，有多少个非零项呢，其实非常非常多的，就是因为中间这个部分它的嗯变换，它是不是特对对应关系还是有点复杂。

如果说呃你在二维的时候，首先你这个u它对和周围的，它其实是和周围的三乘三个呃，节点是有关系对吧，你这个u中间这个和周围的3x3，然后呢，周围33里面每一个节点都有两个component。

所以你这个二维里面就有18个entry，三维里面就更复杂了，虽然三维里面，你这个就是每一个相当于每一个boxo，他会和周围的八个包呃，就就每一个顶点它会和周围的八个box有关系，八个boxo是什么。

八个voxo就有3x3乘三二十七个顶点，你再乘上一个三对吧，因为你每一个顶点上面有三个很ent，所以就81个81个能no zerries啊，你要把这个81个项目都写出来，还是不容易的一个事情。

这k到底长啥样的啊，我们接下来就揭开它的神秘面纱，你会看到他呃，这个k我刚才说它还是非常复杂的，我们就看一个大家一般会讨论的是什么，会讨论一个element，它的这个k e matrix。

ok那么这一个问他这个k matrix，哈哈看起来不是特别好，他首先你得去呃，它有八个不同的项，然后这个八个不同项，关于这个呃排列组合得到以后，得到一个8x8的这么一个矩阵，这个大家一般都是打表的。

我是很久以前手动推过，这个88矩阵里面的一两项已经推了，我嗯挺挺累的，然后一般大家都是呃你抄我，我抄你就把这个88的矩阵啊烂熟于心，也不是烂烂熟于心嘛，就是塞在你的代码里面。

然后呃随时在硬盘里面某个角落能把它找出来，这个地方其实要讲起来就非常非常的复杂了，我这边就不细讲了，但是大家要记得你从u到西格玛中间，这个非常long trivial的这么一个线性变换。

导致了这个8x8的，对于每一个元素的刚性矩阵，有的每个元素钢琴矩阵，然后我们就可以组装出来一个真正的大的，这个k矩阵呃，这个组装过程呢大家可以去，如果你会matlab的话，你可以看看有一个99行呃。

拓扑优化实现用matlab写这么一个文章，然后它是一个非常非常容易上手的一个文章，当然他也会介绍徒步优化，那么既然都讲到拓扑优化，我们接下来就简单介绍一下，拓扑优化是啥，那么我们这边介绍一个。

最最最简单的拓扑优化问题，叫做minimal compliance，比如说呃我们想最小化什么的，最小化它的compliance它的定义是什么，定义是什么，什么叫compliance。

compliance就是呃u transpose乘上k再乘上u，其实就是说某一种形变的能量，或者说这个loss function在如果做gpin的话，你可能更喜欢lost方式。

这个这个称呼他就是mesh呃，你这个物体的形变大小，然后呢我们当然还要满足q等于f对吧，你的这个dispacement得和你的load bor得是对应的，然后你还只有有限的材料。

也就是说你这个空间里面有很多小格子，然后这个小格子里面每一个格子可以是zero，也可以是one，就是可以是空的，也可以是填充的，然后，你就用这个0~1之间连续的数目，来表示这个格到底是有材料还是没材料。

然后大家再用一些trick，让这个值要么是接近零，要么是接近一，把它推到两级上面，这样你就可以把一个离散的优化问题，把它表示成连续的优化问题，当你有个体积限制，你不是说整个空间都得填满。

如果说你要最最小化上面这个东西，不给你几个限制的话，最简单方式是什么，最难就是你把整个空间全部填满，但是你在有体积限制的情况下，你就得想想办法，ok那么我们说了这么多啊。

我们还是看一个demo吧，这个刚才这个例子不是很直观啊。

刚才列了一坨公式，可能大家非常confused，我们来看几个。

看一个拓扑优化的一个演示，找一个白一点的背景。

这，嗯看起来更对，ok好，那么top优化在同一学院里面还是挺常用的，为什么呢，就是因为它好看，倒不是说他做出来的东西呃，当然在工程里面，大家可能更喜欢top优化做出来的呃，这个结构它会省材料的。

在图形学里面很多时候大家就是图它一个好看，但是很多时候由于这个拓扑优化，现在有3d打印的技术，是有3d打印以后，这个突破优化一下子就变得非常的popular，因为这些他优化出来。

材料往往会非常非常的复杂，那么你在以前你用以前的这种增材制造的方式，你往往造不出来，但是你用3d打印以后，这个就非常容易，你直接呃任何材料都可以打印吧对吧，ok那么这个里面我们在干啥呢。

我们这个里面其实是在做一个优化，我们把它的右下角给它固定住，注意右下角固定住，然后呢左下角我们给它呃，或者大家可以脑补一下，这个其实是我用了一个对称的技巧，他这个其实只是结构的右半部分。

因为左半部分和它完全对称了，我们其实就是假设这边有一个点，这边这个点和这个点都是固定的，然后中间给它加一个向下的力，然后然后我就问你，如果就给你这么多材料，你不能把这个空间全部填满。

你只能有空间的50%的volume是呃，可以放材料的，你怎么样去设计你的这个结构，使得它的这个complex，其实也就是说这个函数的值可能尽可能小k，那么呃，祝大家应该也可以在网上下载到这个程序。

你看到这个这个程序是一个交互的程序，你可以调它的分辨率，如果你分辨率高了以后，你近段时间也就会相应的变长，但是你得到的结果也会更好一些，我说的更好，也只是说它看起来更复杂，看起来看起来更好啊。

一般来说有的时候大家会用这个objective，来去作为结果好坏的一个mesure objective，其实就是对应上面的这个lp的这个l ro。

ok那么限于时间关系啊，这个拓扑优化我就不继续讲太多了，大家可以去玩一玩这个demo，这个demo里面有好几种边界条件，刚才我讲讲的是一个比较简单的把脚上固定，这边我们是把呃左上角和左上角和左下角固定。

然后给这边加一个力，这个这个例子好像叫做n b b，并还是叫m m b b，我不太记得了，好像叫m b b b，然后它就是相当于你固定左边两边，然后给你可以想象，如果是左边是一面墙。

然后你要呃优化一个生出来的结构，使得这个伸出来的结构呢啊，最右边这可以挂一个衣服什么的啊，其实还是比较比较直观的，你可以看到它基本上优化出来，这么个像钳子一样的结构，还有这个边界条件三。

边界条件三是一个桥梁结构，这可能分辨率有点高，一下子算不出来，ok这个是干啥的，这个是把左右左下角，右下角给它固定住，然后给上面每一个点都施加一个向下的力，实际上就是一个桥梁的一个设计的一个东西。

但是我不知道现实中有没有人真的这么设计，桥梁，可能会有有些建筑师他会喜欢用徒步优化呃，来设计这些建筑，其实有真的，有些建筑是头部优化设计出来，大家不要觉得这个头部优化是一个，理论上的东西，很多。

由于它设计出来这个结构往往非常具有美感呃，很有这种叫做生物质感，所以很多设计师最近也比较，最近不是很流行这个环保仿生这些设计，所以很多设计师啊，他会用徒步优化来去设计这个建筑，这个城市有的时候啊。

sob有限元还是要费一些功夫的，你可以看到他这边在非常卖力的去迭代，看到他这个由于我们用的是ctrl redis，没有加预条件，所以你可以看到它收敛速度并不是很快，特别是当你的这个结构不是很清晰的时候。

大家如果要看这个minimum compliance，拓扑优化其实只有两个关键词，一个叫做sip呃，一个叫做oc，你只要记住这两个关键词，然后到那个文档里面去找，这simple和oc是啥就可以了。

这个其实算法和之前推的有限元的比起来呃，就和非常容易理解了。

它没有那么多非常复杂的东西，它其实就是一个很简单的一个优化问题，ok那么说到这个头部优化，也不得不提，我们之前做过一个非常有意思的工作，就是说呃，如果说我们呃把这个图优化的分辨率，推到非常高会怎么样。

我们刚才演示的top优化，它分辨率其实非常低的，他只有大概100x200，这个分辨率图形学里面大家很少用这么高，这么这么低的分辨率，大家一般都喜欢，至少是几百乘几百啊，至少是大几百。

或者是这个几千是几千人这么个分辨率对吧，那么特殊优化，由于它是要算一个f e m，所以它的分辨率之前一直都很难做的很高，因为你要solve fm还是挺费时间的，你会看到之前，首先他这个stl非常大对吧。

和普通比起来，然后非常复杂，然后，呃他还你要给他设计一个presugitioner，还不是特别容易呃，ok然后所以说我们就做了这么个工作，这个工作其实它其实目的非常简单。

就是说如果说我把这个分辨率调到非常高，然后只给你一台电脑，只给一台电脑，我们去做这个拓扑优化啊，解fm，然后我们能得达到一个什么样的效果，这个视频我传到哔哩哔哩上面，大家可以去看一看。

应该是刚审核通过的，昨天晚上传的，它其实里面核心部分就是一个高性能的f m server，然后当时我们找了一台机器，找了一台有512g内存的机器，这样我们可以在一台机器上面，就用cpu去做有限元求解。

然后推到一个很高的分辨率，这个其实有应该是有超过10亿个，是不是11，这个是呃，对应该是10亿对10亿个呃，超过11个这个boxes，然后你会可以想象这个里面。

你如果用光把这个k matrix给他存下来，就不够存，那怎么办，你得做matrix，也就是说高精度的情况下，大家往往会考虑matrix的搞法，因为你没法把这个矩阵存下来。

你还得设计一个很好的precondition的，要不然你这个收敛其实就收敛不了，然后呃，你或者你要很多很多次cg iteration才能收敛，就效率非常低，呃，这个其实根本质上就是一篇高性能计算的文章。

有我看到有同学呃说到这一点，他其实呃和graphics相比，和传统的graphics simulation比起来，它更像是simulation里面加高性能计算，ok这个是上面显示的。

这是一个鸟嘴结构的拓扑优化，这个是鸟嘴结构是干啥呢，就是说我如果把这个面固定，然后我给他这个鸟嘴的某些部分家里，因为你这个鸟嘴得是一个有用嘛，你鸟嘴得够硬才能去吃一些果子呀，什么啊。

鸟一般是吃什么坚果对吧，鸟嘴里够硬才能咬得动这个啊，花生啊什么东西的，所以这个设计上面，你是要设计一些有一个条件的，ok你就可以看到它这个其实得到一个仿生设计，其实你可以去看一看呃，有些鸟它的化石。

它其实内部结构可能和这个还是有点像。

不能不能说完全像，毕竟这个是在计算机里面进化出来，它不是在自然条件底下进啊，进化出来的呃，ok那么其实今天我们就讲的内容就差不多了，我来看一看有些问题哦，有同学说拓扑优化的那两个关键字对。

其实这个拓扑优化它其实关键字就是就是这样，你只要在这个文我的作业包里面有一个pdf，但这个作业大家只是有空的话就可以做一做，没空的话就不用做，那个代码是c加加写的。

因为这个是早期的时候太极还是一个c加加库，然后我后来实在没有找到时间，去用新版的python的太极重写了，下次一定这次真的是没有时间了，然后他是一个c加加的作业，然后我们在m i t交通信学的时候。

也会用到那个作品，当时是你会看到它里面有一个2万行的图文件，然后有一个200行的拓扑优化的一个程序，那个代码可以在各种平台上，不用装虚拟机就可以运行，然后那个文档里面那个包里面有个pdf。

那个pdf里面会介绍sip o c，然后还有手把手的教你，怎么把这个代码从一个不word的这个程序，把它呃，优化到一个逐渐实现到一个能word的头部优化，最后你能看到就是我刚才展示的这个。

这个程序的效果，有同学说f m能不能结合multigrade，对你如果要看fm加mt grade，那你一定要看一下这篇这篇文章，这篇文章是fm加mtv的一个，我觉得还是不错的一个文章。

当然这个mt贵不是像普通的mtv那么容易啊，这个mt贵的还是非常头疼的，要把它写对真的很不容易，你可以看这个代码其实是开源的，但是一个c加加的一个代码呃，为什么不用gpu，因为gpu没有这么多的内存。

gpu它的memory可能只有现在可能只有22g，24g，但是我们整个系这个这个信号跑起来，大概要用到300多g，或者将近将近400g的内存，但你不要觉得300多g400 g好像很多。

其实之前有一篇nature paper，也是做这个高精度的table depremization，然后他用了我记得用了500台电脑，还有什么8000个cpu core才把它算出来，我们只用了一台电脑。

50个cpu块就把它算出来了，分辨率和它也差不多啊。

GAMES201：高级物理引擎实战指南2020 - P7：Lecture 7 混合欧拉-拉格朗日视角（1） - GAMES-Webinar - BV1ZK411H7Hc

好那我们正式开始嗯，今天我们来讲一讲混合欧拉拉冈人的视角，而这个是一个听起来很复杂，其实非常非常简单的话题，我们这一讲啊，其实没有任何公式啊，这样非常好，所以我直接都用keynote做slides。

没有用雷tag呃，下一讲我们会讲一些讲一讲这个背后的理论，这一讲我们就专门讲他的呃高层的一些想法。

今天呢我们主要讲一讲，整个欧拉拉格人视角的一个overview，然后啊我们具体介绍一些例子，遗留格式，比如说particle excel，fm。

particle excel和polynomial particle excel，呃，他们分别简称pick a pic，polypc，这三个pig方法是中医学里面非常非常常用的，混合欧拉拉格尔方法。

然后还有一类方法呢叫做flip，flip，其实也非常非常常用，今天我们呃除了讲这几个呢，我们最后还会讲一讲m p m，然后下一讲呢我们来讲m p m的理论和实践。

以及m p m里面的这个constitute model，好像中文叫本构模型，然后我们会讲pm里面的拉格朗日历，以及影视，m p m时间积分和一些高级的采集feature，这是下一周的内容。

这周呢我们就讲一讲一些比较基本的内容，好我们先来回顾一下什么是拉克人视角，什么是欧拉视角，他们其实是连续介质的两种视角对吧，那在拉杆人视角里面呢，你想想呃，你可以把自己想象成一个随随波逐流的小船，对吧。

呃你每时每刻都在问自己，我的位置和速度是啥，你是随着你模拟的材料一起移动，这个就是拉格尔视角，一般来说在拉杆视角里面，大家喜欢用party，但并不是拉格尔视角里面就一定就是party go。

还记得我们之前的有很多有限元的模拟，它也是拉格尔视角对吧，但是它就是用的是一个flash，ok那么欧拉视角呢，他和拉格尔视角是不一样的，欧拉视角里面的这些所谓的传感器，或者你自己的观察的位置。

它是从来不移动的，然后你每时每刻都在问自己另外一个问题，这个问题是对于我自己这个位置经过的材料，它的速度是啥对吧，那这个就是拉格尔视角，其实还是非常容易理解的，你可以认为它就是一些插在水里面的桩子。

如果你在模拟这个水的话，这个桩子是从来不移动的，这个克拉管视角是形成了鲜明的对比，还记得那句话吗，呃拉港人视角是随波逐流，在统计学里面，大家两个都经常使用纯拉格朗视角。

和纯呃欧拉视角都是非常常用的模拟方法，在拉格朗的视角里面，大家一般来说会用particle，然后在欧拉视角里面，大家一般会用一个grade，或者用一个mesh，那拉冈伦视角里面也可以用mesh。

也可以用，但是相对来说比较少一些，拉格尔视角里面典型的一个呃，一系列方法是基于粒子的方法，比如说s p h啊，然后position based dynamics，里面很多东西都是用八个人视角来做的。

比如说你看左边这个是nvidia做了一个工作，他们所有东西全部都是用例子来表示的，你可以看到拉格伦视角它的很大的一个好处，就是它用粒子来表示以后，它可以表示出各种各样的呃形状。

然后各种各样的材料也可以用这个来表示，欧拉视角呢，这边右边是一个很经典的一篇烟雾模拟论文，也快有20年了，他是用一个背景网格来表示，这个厂里面的速度和density，这个烟的浓度之类的。

各种物理属性都是在纯grade上面去表示，没有粒子，那么很多同学自然要问，到底拉格朗视角和欧拉视角哪一个更好呢，那么一般说到这种问题啊，我们都得先定一个标准对吧。

嗯那有好几个key facturers是关键因素，我们得去考虑，那其实我们在这边就是在呃define，什么叫better，better有好多种层面对吧，如果一个人跟你说，这个算法每个层面都好。

那你可能得想一想他是不是在忽悠你，因为很少有这样的算法，一般算法也就是在几个层面比较好，你要根据具体的情况来选择使用算法，真正的在每一个层面都比较好的算法，是非常非常少见的。

ok那么我们这边要考虑哪一个更好的时候，要考虑哪些因素呢，首先我们要考虑的是物理量的守恒对吧，你这个模拟的方法你至少得实现动量守恒，然后在要求高一点，可能得角动量守恒对吧。

那你还可能得考虑你的体积是不是守恒，如果你在模拟的是incompressible的东西的话，比如说你在模拟嗯不可压缩的流体，那你希望它的体积能够守恒，然后你还希望它的能量能守恒。

比如说你这个机械能不守恒，然后你别看这几个守恒量看起来非常tribute，但其实你真正的要把它做到的话，还是非常不容易的，你很可能设计一个scheme，自己弄着弄着，你发现诶好像角动量不守恒了。

然后呃一般来说减重量不守恒，很少来比较少，是爆炸的情况，一般来说是慢慢慢慢就没了啊，这个其实和能量的耗散其实也非常相关，那接下来呢我们就要考虑它的性能性能，在现代计算机体系结构上面的性能。

基本上就是两个factor，一个是parallem，一个是quality，你这个方法容易不容易并行，然后这是计算上面呃，访存上面是你的这个算法，在访问内存的时候。

是不是有比较好的special和这个temporal locality，它的时间和空间的局部性是不是好，因为在当前的计算机体系结构上面，你的flop就是你的这个浮点数计算，是一个非常廉价的东西。

真正昂贵的，不管是从时钟周期还是从能量上来讲，都是你从主内存移动到寄存器，这个移动内存或者叫保存的过程，所以考虑paralysm和locality，你得嗯和当前计算机体系结构进行结合。

然后来考虑这样才能达到一个比较高的性能，最后呢你可能还得考虑考虑计算实现的复杂性，对吧，你这个算法有可能各方面都好，但是就是实现起来太复杂了，没人能正确实现出来，你别笑啊，其实有很多算法都是这样。

他好像呃看起来非常简单，你真正要实现起来还是非常非常难的，不是所有的算法都能被你，99行或者88行给实现是吧，当你这个复杂程度高了以后，很有可能复杂的程度高以后，他就不守恒了，或者性格就变低了。

为什么呢，为什么不守恒了，也没人能实现对，那当然就不守恒了，为什么性能变低了，因为你算法变复杂以后，你要去优化它就非常非常难啊，你随便你找一个简单的算法，你教给我，我给你重写一遍，快十倍。

一般都不成问题，但是你如果给我一个非常非常复杂的算法，那我就要考虑一下我要不要帮你做这个事情，因为这工作量可能成本太高了对吧，我就呃即使一个非常擅长优化的人，他可能也不愿意花这个时间去帮你优化。

因为太复杂了，ok那么讲完了这么多呀，我们定义好的一个标准，接下来我们来看一看，我们的混合欧拉拉格尔方法，为什么我们要用混合欧拉拉梗式方法呢，当然是因为我们会让他们希望，让它变得better对吧。

让它变得更好，那我们想想看，在一个fluid solver里面，通常我们有两个component，对大家还记得吗，如果是一个incompressible fluid solver的话，一般就两部对吧。

一个是the action，一个是projection，action，就是说你要去移动你的流体的这个场，projection是什么呢，projection是你的速度长，如果说有散度分量的话。

我希望把这个散度分量给他投影掉，也就是说我希望去enforce他的incompressibility，我希望去满足它的不可压缩性，大家在之前的实现中应该逐渐也发现呃。

action和projection系呃，喜欢的表示其实是不一样的，欧拉网格呢就特别擅长projection，为什么你在欧拉网格上面还记得吗，在欧拉网格上面，如果我们用一个均匀网格的话。

那很容易district，还记得我们的五点播送stl吗，对吧，就五点泊松的模板其实就是中间是二维，里面是中间是-4，上下左右各是一，或者你取一个负号，反过来也全部都没有关系。

均匀网格上面聊起离散拉普拉斯算子，是非常非常容易的，同时呢你如果对于一个网格上的节点，你要去查找他的邻居也非常容易对吧，那你这个在内存里面都是均匀排布好的，你只要呃比如说你如果是i g的话。

那你上面的就是i j减一下，面就是i j加一对吧，这个其实是非常容易实现，neighborhood，look up，非常容易实现，然后呢你如果要去做preconditioner。

由于你的结构是一个均匀的网格，这边你要做precondition也非常非常容易，很容易你就可以上一个geometric multic grade，这个几何多数网格法去做preconditioning。

这样你收敛也会非常快，所以你会看到欧拉网格它其实在做projection的时候，它有非常非常多的优势，但是它有一个很大的缺点是什么，他在做the advection的时候很容易出问题。

还记得我们之前呃讲了好几种vection的方法对吧，不管你是用sama runner也好，不管你是用呃micromic或者b p c c也好，呃你可能有一些方法比另外一些方法是好一点。

但是他总改变不了它的数值耗散的问题，可以看到我们右边这两个图，假设我们一开始的时候，这个太极图标是这样，当他转了几圈以后，自己往往就会变得非常的模糊，他会丢掉能量和几何上面的细节。

这个就在模拟的时候会非常的不尽人意对吧，因为你模拟的模拟人，你这个能量就跑了，跑没了，那这个就呃在流体魔力里面，往往这个导致的结果就是说你的呃，流体看起来非常非常的黏。

或者说一般大家叫做two biscu，它的联性太高了，然后同时呢你可能有一只，一开始是一个兔子形状，你这个模拟就模拟它就变得越来越糊，然后他就不知道变成一个什么形状，最后可能就变成一团球了。

那这个就呃看起来效果就不是很好，当然你如果说你要考虑准确性的话，那准确性也不好，ok那么刚才讲的是欧拉网格，接下来呢我们就讲一讲拉格朗的例子，那么粒子它非常非常擅长一流对吧，我要移动一个例子。

或者只要干啥，我只要把他的位置加上dt乘上速度，或者用一些更高级的格式，我就可以简简单单把这个物理量在空间中移动，同时呢他或许还会非常呃conservative，比如说你如果把质量track。

用例子去track你的质量，那你的例子由于你从从头到尾它都不生不灭，所以你的质量自动就是守恒的，动量呢你也很容易做到守恒对吧，那如果你用规的话，那你移动以后，你可能得做一些额外操作，才能保证动量守恒。

这个是它的好处，它非常非常擅长，但是你如果在例子上面做projection，那你这个就非常头疼了，为什么呢，因为你粒子在空间中的分布啊，你的呃它是不规则，他可能你中间有一个例子。

那你怎么定义他左边的例子呢，它左边可能有好几个例子，你很难定义清楚，他要离散化，它就非常的tricky，之前哦，我们很多同学实现s p h里面这个shape function啊。

这个叫kernel function，它结果好多好多种变种，为了在粒子上面能够实现离散化，这个实现起来还是有点头疼，同时呢你要去找你的邻居，你要去apply standal的话，你得找你的邻居对吧。

那你的查找邻居这个过程也非常非常的复杂呃，一般来说大家都要用一些比较复杂的数据结构，比如说一些hash reen，或者一些呃，各种各样更fancy的数据结构来实现邻居的查找。

一般在纯particle的方法里面，邻居的查找以及在邻居上的操作就是它的瓶颈，所以例子它并不是真正擅长，projection的一种表示方式，那么问题来了。

既然alerin green欧拉网格擅长projection，例子擅长the action，我们有没有一种方法能够把两个表示各取其长，然后又能避免其短，把他们smart combine。

就我们把他们聪明的结合在一起呢，答案是肯定的，这个就牵扯到我们今天的呃，混合欧拉拉格式方法，既然我们要做到科举学长，那我们得两个表示都得有对吧，那么一般来说呢，在混合欧拉拉格人方法里面。

他的欧拉grade和拉格朗日particle，他们并不是同等地位的，一般来说大家是呃，把拉格尔的例子作为一等公民，欧拉网格作为二等公民，为什么呢，因为你大部分信息啊，主要还是存在拉等人的例子上面。

the grid呢一般只是一个附属的，用来计算中间结果的东西，ok，那么明确了拉格的例子和欧拉网格的地位以后，我们呃就可以设计整个算法，首先呢我们有一些拉杆的例子对吧，它存所有的状态。

然后我们要一般来说在这一切方法里面，它的套路都是第一步，我们先做particle to retransfer，我把网呃例子上面的数据传输到网格上。

这一般大家叫做p t g particle to grid，简称p g，那么做完p to g以后啊，你的网格上面就有一些物理的信息对吧，那我们在网格上面就可以做一些grade的operation。

网格操作，这些great operation有一般有啥呢，一般有呃pressure projection，你希望去呃把你的呃，速度的散度分量给他投影掉，去enforce你的速度的无伞或者不可压缩性。

也有可能在grade上面去做一下boundary condition，但boundary condition也可以在party上面做啊。

我们这边只是说grade上面的boundary condition，边界条件，然后gre上面一坨操作给他做完了，然后我们就可以做great to particle transfer。

我们再把信息在gre上面走了一圈以后，我们再把它重新拿回来，然后这样我们就又把信息到了particle上面，最后呢我们在particle上面做一些更新，包括移动我们的particle。

然后或者呃更新材料的一些constitute model之类的，比如说更新deformation，gredient，更新它的形变梯度以及更新它的体积啊，一些各种其他的属性。

这个其实就是整个混合欧拉拉格尔方法的呃，这个套路其实还是比较清晰的，他就是有既有欧拉表示，有拉格朗表示，我们用transfer方法，把信息在两个表示之间来回倒，那么这个来回倒的过程。

它其实也会引入一系列的问题，包括你来回倒会导致信息的损失，来回倒会导致计算量变大之类的，各种各样的问题，我们在后面都会逐一解决。

那么我们现在就来考虑一个最最最简单的呃，混合欧拉公式方法叫做particle excel，一般大家简称pc或者读作pick，那然后这个party insult方法，其实就是一个非常非常古老的一个。

混合欧拉拉控制方法嗯，在pk里面，他的p图记的过程其实也非常非常直观，就是对于每一个particle呀，我把他的各种物理信息给他，传输到它相邻的若干个格点上面，这边我用一个二维的二维的例子。

然后他的呃传输范围是3x3的great node，然后传输的量是什么呢，一般来说最常见的就是速度，偶尔呢大家也会传输温度或者力，m p m里面就会传输力，然后这个是particle degrade。

然后你在gree上面做完一些操作，你得做great to particle对吧，那你其实每一个例子，从周围的3x3这么个网格上面去收集信息，然后收集的一般来说就是呃，速度或者温度之类的物理量。

那么在传输的时候啊，一般大家会遵循一个原则，就是说这个例子并不是均匀对待它的邻居的，或者说呃是不平等的对待他的邻居，当一个邻居节点对他更近的时候，他就会有更多的importance。

这个这个大家应该也能理解，就是说呃距离更近，它的重要性就越高，那么一般来说呢这个importance是怎么定义的呢，大家会用呃一些和函数，或者叫做你把它叫做差值函数，也可以这边叫差值函数。

可能有点奇怪啊，你可以叫它，这个就可以叫它是一个kernel function或者核函数，然后这个核函数呢一般来说大家的搞法是，我在x方向上有一个和函数，在y方向上有一个和函数。

然后这两个核函数呢把它乘积起来以后，就得到在呃每一个格点上它的权重对吧，和函数一般是怎么选择呢，呃和函数其实混合欧拉拉公式方法里面有很多，很多种选择的方式，然后我们今天介绍三种最最最简单的核函数啊。

这个叫b splenko，然后一般来说呢，大家呃最常见的俾斯麦kernel是什么呢，其实你如果做线性的这个差值，就或者这个帐篷函数，它其实就是一个最简单的linear belan。

然后呀你如果希望它长得更smooth一些，你可以用quadratic belin，quadratic belin就是下面这个图里面的红色的线，零点就是这个黄色的线，然后你可以用cube cubic。

它就是下面这个蓝色的线，你可以看到，当这个b splay它的变得复杂程度变高以后，它的计算量其实也稍微变得更高一些，如果你是零六的话，你就要取一个f a b s对吧，你就要取一个浮点数绝对值就可以了。

但quadratic呢你就得算平方操作，cubic就有三次方的操作，然后一般来说大家用linear和quadratic，用的比较多，cubic实际上用的相对来说比较少一些。

因为quetic一般已经足够smooth了，linear有什么问题，linear它不够smooth，然后有的时候会导致你的模拟系统，就不是特别稳定，所以project一般是大家用的最多的。

特别是在m p m里面，在flip里面大家在流体里面大家会用liner，然后后面我们会讲，ok那么我们来看一看particle excel，particle degrade怎么个写法啊。

其实也非常简单了，这边假设我们是用一个sales center grid，也就是说我们的呃，节点呢是放在每一个cell的中心，那这是什么意思，也就是说我的每一个cell，他的顶点的位置是dx的倍数对吧。

比如说这个点就是零零，然后这个点呢是dx 0，这个点呢就是2d x0 ，那么所有的呃，所有的这个node它的位置其实就是必须是多少，加上0。5的dx，比如说这一个note，它就是2。5d x0。

5dx横坐标2。5呃，纵坐标0。5，ok啊，然后如果说我们用这样一个规则规则的情况，我们对于每一个例子啊，对每一个particle in x，这个x就是所有的particle的位置。

那么我们得计算它能够它的3x3的范围里面，左下角的那么一个节点坐标是什么，对吧，那么怎么来算呢，首先我们把节点的位置，我们so sorry，我们把我们把粒子的位置乘上一个inverse dx。

也就是说我们把它除掉dx，那么我们其实就是把物理空间里面的连续的位，置，把它转化到了网格空间里面的呃，这呃这个离散的位置对吧，我们除掉dx粒子，除掉dx，然后减0。5，为什么减0。5呢。

注意我们这边格点的位置是和呃，它的位置一定是0。5，dx他一定一定有这个0。5dx分量在里面，所以我们把这边减到0。5，ok然后呢我们算这个quadratic display的这个权重。

然后这个权重算出来以后，注意我们这边其实是一个所有的这些操作啊，全部都是二维的向量的操作，所以我这边其实算了多少多少个值，我算了六个值对吧，三个值是x分量，三个值是y分量。

然后我们把这个权重x和y方向上的权重，分别算出来以后呀，我们每一举它的范围里面的3x3的node对吧，node就是格点了，然后我们对于每一个note。

每个note它有一个相对于左下角的node的offset对吧，这offset其实就是什么，其实就是ig对吧，我们i j是从0_{3呃，呃i从0}3，这个从0~3权重是啥呢，然后i是x方向上的偏移对吧。

那么我们取权重的x分量，然后第一个vector，然后取权重的d y分量比这个vector，然后我们就把x分量上的权重，和y分项的权重给它乘起来，得到了它这个节点的权重，然后呢我们就得把信息。

把这个particle上面的信息给他写到节点上面，注意我们这边每一个节点它存了两种信息，第一种信息是你真的物理量，这个真的物理量就是你的权重乘上你的物理量，给它叠叠加上去，那另外一个信息。

另外一类信息非常非常重要，是什么，就是你的权重只有你的权重，为什么呢，因为你有些地方可能例子多，有的地方可能粒子少，所以呢你就希望去呃把我的权重也算一算，然后这样我们就可以去做。

用权重呢去做一下normalization，这样实现一个守恒啊，这边可能稍微有点需要理解一下，其实就是每一个从每一个节点的视角来看，它的物理量是什么，它的物理量就是周围的particle。

立particle的物理量的加权平均，那你加权平均，你加了权，那你总得把全的和算出来，然后把全除掉吧，总总有一个全的normalization的过程，所以呢他也需要去算一算权重的和，除了物理量的核以外。

ok那么刚才说了，我们这个有一个权重的和，那么其实就是gram这个物理量对吧，然后我们在做完particle degrade以后，我们得做一下great normalization，对吧对吧。

物理量硅上面物质要除以这个权重，首先呢我们得我们只需要处理，权重大于零的grade就可以了，如果说它的权重是零的，这个根本就不在particle能够触及的范围之内，我们也就不需要去处理它。

那么对于所有权大于零的例子呃，这个网格点呢我们就把它的速度除掉，它的权重对吧，这个速度其实是有加权平均，所以我们要除掉甲醛和这个非常容易理解，那么做完那往往你做完normalization以后。

你得去做各种像pressure projection啊，之类的各种操作对吧，那么我们这边就呃忽略掉pressure projection，或者忽略掉我们节点上面受力的计算了，我们直接考虑呃。

grade particle，我们直接考虑它的advation的部分，和刚才其实规律particle刚才其实一模一样对吧，我们还是要算一下它左下角节点的坐标，然后算它的权重，这个和刚才一模一样。

然后我们局部啊有一个叠加变量是这个v，就是accumulate呃，我们把所有的网格上面的v都加圈以后，叠加到局部的这个v的一个和对吧，其实非常容易，这个其实非常容易理解，其实就是每一个例子。

把周围的节点上面的速度给他收集一下，然后重新收集到这个例子上面，那么算出来这个新的v啊，我们就可以呃，我们可以把这个x用我们先移动这个例子，我们把x用上一步时间的v，然后乘上dt该累加上。

当然你这边也用新的v也是可以的，这边其实是一个呃显示欧拉法，那如果用新的v呢，然后我们用新的b去更新粒子的b，这应该还是总体来说是非常直观的，那么刚才我们讲的其实只是particle excel。

那么如果说这个纯particle excel其实没什么用对吧，因为你只是在移动例子，你没有做什么有意义的操作，那么particle excel呢，往往可以和基于网格的泊松求解器联合起来。

然后得到一个流体模拟器，ok那么它的这个simulation cycle是长啥样的，其实和刚才非常简单，我们在particle，particle，particle to grad，操作里面。

去把速度从particle transfer到grade上面，一般这个这一步大家会用scatter，然后great normalization，然后呢我们在规则规则上面去解它的pressure。

解完pressure以后呢，我们把这个呃pressure apply在你的速度场上面，apply到速度场以后，我们用g9 p的操作，从你的速度场去gather information。

然后从网格到粒子上面，那么然后接下来我们可以根据呃，这个无伞的速度场去移动一个例子，注意它这边已经变成了，由于我们有一个pressure projection的过程，这边新的速度上就是无伞的五散。

主要是为了什么，主要是保持体积对吧，你是incompressible fluid，这个时候你移动例子的时候，你可以用2k two three four，你都可以愿意用一个高阶格式的话，会更准确一些。

那么我其实有一个demo，这个是我老早写的，ok我们来看一看这个demo，可以稍微用命一下，啊这个是我老早写了本科手写了，换了1080ti的机器，现在我们可以随随便便就上最高分辨率了。

我们来稍微调一下这个体积，subs steps给它调大哦，我还不能缩小，这边还有个重要的参数，我得调整零，这个其实就是一个particle excel simulation。

注意你做particle excel的时候，你要记得把这个demo里面的flip blending，给它调成零，因为如果你的这个flip blending不是零的话，它会把一些flip的成分也包含进去。

这边我后面会解释什么是flip，ok你可以看到它体积不是很守恒，是因为我们的，这个sober他的迭代次数不够多，这个地方没有用，country ingredients。

这个地方由于为了gpu上面实现方便，我们是用的jo be对吧，大家可以在自己机器上玩一玩，你可以看到这是一个party excel simulation，它是流体的，它总体来说看起来不是特别好。

它看起来会稍微有一点点，如果你特别是比如说，我们把分辨率稍微调低一点，这样它可以，你可以看到它很多的这个涡量啊，它这个旋转这个旋涡转着转着模拟模拟就没了，特别是你分辨率越低的话，它这个数值耗散就会越大。

比如说这个128x128，它看起来非常灵，就像一团糖浆，或者说呃也不能说蜂蜜也分，当然没有蜂蜜这么甜，可能是蜂蜜和水的混合物那种感觉，一般来说大家看到这种情况就会说它过于一点，然后不太会使用。

所以这个party excel有一个什么问题，他这个能量耗散就比较严重，一般大家很少直接使用party incel，那么为什么会有能量耗散呢，我们来看一个最最最简单的。

particle incel的一个demo，这个大家代代码包里面其实是有的，叫pick versus a pig，我们可以来跑一跑，我把它启动一下。

ok这是一个最简单的party political excel的demo，然后我们现在是用的，你看底下有一个scheme，你可以用a来切换这个scheme，那个按a的话，我们还没有讲到apex。

所以我们现在就只考虑这个pk，如果说我们把这个速度场，初始成一个translational的速度，成比如说向右的速度上啊，你会看到它移动的是非常顺畅的，注意啊。

我为了防止它跑出我的simulation demin，我就再把它的速度在边界上面做了一下clamping，把他位置climb clip一下，防止它跑出去，你看到他的如果是平移的话，它是非常顺畅。

但是呢如果说我按一下r，把它的初始速度设置成一个rotation的话，你会你会发现它转着转着它就不转了，这个就非常头疼，因为你在模拟流体里面，你很多时候有这个窝对吧。

他war本质上就是流体在做主tation，那他不转了，你这个就很影响流体的效果，那类似的呢，我还可以把它按，我还可以按键盘上的d把它设置成dation，设置，带雷神以后，它就是一个拉长的效果。

你会发现啊，呃这个dation它也有问题，他这个拉着拉着它就自己就不移动，理想状况下呢，它应该是我可以看大家看一下，理想中的话，这个dation应该是不停的在拉伸自己。

ok那我再切换回我们的这个呃particle in ca skin，那么还有一个sharing这个运动，sharing就是这个剪切的运动，这个剪切的运动它其实也不能保持，他就模拟着模拟着能量就耗散掉了。

那么为什么会有这种现象呢，这个就是你得考虑它的，great to particle的这个传输的时候啊，它就会有这个信息的丢失，为什么呢，假设我们只有一个例子，kernel support是3x3的话。

那么我们有33x3=9个格点，每一个格点上假设我们传输了，我们只传输速度啊，那只传输速度的话，每一个格点上有x分量和y分量的速度对吧，那它其实有18个自由度在格点上面。

但是你的particle上面有几个自由度，你八六上面的速度，你只有x方向，x方向上速度和y方向上速度，你八六上面只有两个自由度，那你这个传输的时候，从18个自由度一下子变成只有两个自由度。

你是不是就得丢失信息对吧，那么怎么样来解决这个问题呢，有两套解决方案，第一套解决方案呢是，我们不要只在例子上面保留两个自由度，我们可以传输一些信息对吧，我们粒子上我们可以在网格上面看一看。

他是不是还有旋转的这个分量啊，它有个角速度，有没有一些sharon的分量，有没有dlational的这个分量，我们都把它记下来，两个方法，这两个都是比较接近的方法，比较呃最近提出来的方法。

epic是2016年seagraph，然后polypg是我没记错的话，2017年的c graph asia，在呃，都是在通信学里面提出来了两套模拟的方法，另外一个呢另外一套搞法呢。

就是说我不transfer更多的信息，但是呢我呃不把所有的信息都通过grade，来传传输，我只通过gree传输一个dota，传输一个dota，这样我们就能在party上面保持一些信息。

ok那么我后面会介绍flip这个方法，那么apex其实直观上来说是非常容易理解的啊，刚才我们看到如果我们只用pc的话，它只有这个constant，它只有这个常量速度场的两个分量对吧。

这是一个v0 和v1 ，我们这边用零表示x分量，一表示y分量，就是呃x轴和y轴，那么实际上我们刚才看到，我们刚才还有很多常见的速度的呃模式，比如说我们有这个c00 ，c0 是什么，是dation。

是x分向上的dational的这个分量，就说它沿着x它的速度不是一样的，你的x越大，它的速度越大，它其实就是在拉伸你的材料对吧，这个叫c00 ，还有c10 ，c10 呢是一个sharing的这个动运动。

它随着你的x的增加，你的y方向的速度会增加，它是一个sharing的一个mod，一般来说大家会叫这些速度场的这个pa叫做mod，然后pk的话其实只有两个mod，apex，它的mode就会更多一些。

ok那么然后还有什么，还有c01 ，c01 什么也是一个sharing，随着y变大，你x方向这个向右的速度也会变大，类似的还有c1 ，c一是y方向上y方向上的dation。

比较复杂，实现起来却异常简单的一个搞法啊，这个他的直观上理解也非常非常容易，我记得我是第一次看到这个apex的论文的时候，我觉得好像非常复杂，另一个时候我大概大三上学期的时候，我记得非常清楚。

因为算法的复杂，这个数学上的复杂程度和实现的复杂程度，反差过大，所以这个算法给我留下了很深刻的印象，呃，大家可以在哔哩哔哩上看一看这个epic的视频，这个强烈推荐这个整个视频啊。

包括这个论文写的都是极其的通俗易懂，而且呃写的可以说这个paper writing以及这个video，它的制作质量都非常非常的高，我强烈建议大家去看一看，你会觉得收获颇丰，即使你不能理解其中的数学。

你光从他的high level idea上面去理解一下也非常好。

之前我们刚才看了这个旋转的这个mod，用pc的话，它就叫流量都没了对吧，但是apex它的particle，degrade和grade particle都是角动量守恒的，然后这个推导比较复杂。

大家可以去看两篇文章，一篇文章是发在西瓜上面，这个推导，大家有兴趣可以看看，没兴趣的话，你不看也没有关系，大家可能看到刚才很多，可能很多同学在这边会被刚才那个推倒吓到了，说这玩意儿实现起来效果这么好。

然后数学上公式有那么复杂的，实现起来是不是也很复杂，其实完全不是，它实现起来只要三行代码就可以了哈，这个呃可能都不是三行完整代码，就2。5行代码就非常容易呃，我们来看这个p two g和g9 p操作。

那p9 g操作呢，其实和刚才的pk几乎是一模一样对吧，那唯一的一个变化呢就是我这边有一个i fine，这个所谓fm particle excel吗，它有一个f一的速度场。

这个fm速度场我们再从particle，转移到规则的时候呃，除了我们去呃转移particle本身的速度，这个平移速度这个v我们还得转移呢，它的一个局部的这个fi速度，i fi速度场是什么呢。

呃fm入场就是你的这个c乘上d p o s，d p o s是这个节点的位置，减去particle粒子得到一个数值，然后这个是一个矩阵乘法，这个具体的数学上推导，我们会在下一讲里面更深入的讲解。

这一讲呢我们只讲具体实现，因为你如果能实现出来，那你就有更多的兴趣去研究它背后的数学公式，一开始就讲数学公式可能比较枯燥，对于我这种学计算机的人来说，当然你如果是数学公理非常好，你如果是学应用数学。

那你可能会先想想推公式，但我相信大部分同学还是比较喜欢，先把它实现出来，再去讲其中的道理，ok然后这个是p9 g g two p呢g6 p的时候，你得去算一算你的c对吧，那算这个c其实也非常容易。

就算c其实你只要去对于每个节点，你把它对c的贡献给它叠加起来就可以了，那这个地方就有一个一系列系数，然后一个有一个这个是output a，这个是两个向量的外籍，这个具体为什么是这样啊，我们先留一个悬念。

下下一讲再讲，这一讲大家就先接受它，就是这样，ok那么可以看到这个公式，可以这个实现可以说是非常非常简单了，那么我们来跑一跑，刚才我们跑的是pick，我们同样的程序，什么样的效果，还记得吗。

pick我不管给他什么样运用，除了translation pick是好了，这个是还是用pig translation pick是好的，我看看dation pick不行，然后sharing pk也不行。

rotation呢pk更不行了啊，rotation可以认为是dation和sharing，组合起来得到的东西，那么现在我按一下a，你可以看到我这边有一个scheme，你可以用a按a来切换scheme。

我们来按一下a，我们现在把它切换到了epic，那这个时候我们看看translation和pc一样，pick translation也没啥问题对吧，我们来看rotation。

你可以看到他这个rotate起来就非常的，相当的流畅，rotation不会出现pick刚才亮的呃，转着转着就能量耗散掉的情况，我们再看一下pick，我们按一下啊，sorry sorry，我们换一下a。

把它换成pick，你看他就不行了，那sharing呢，sharing也没有，完全没有任何问题，打雷神呢，打雷神也没有任何问题，呃呃我们和这个pick dation比较一下。

pk雷神和能量耗散非常非常之快对吧，因为pick它倾向于去smooth你的输入场，他会不断的把你的速度场去平均化，但是dation呢他们思路上每一点都是不一样的，那你在嗯做pig的时候。

你的速度上就会被逐渐的抹平，抹平了以后它就不动了，好可以看到这个epic是像可以说是效果惊人啊，很多时候大家还不满足，那这个时候就会有一个polypc，polpc是什么个搞法呢。

你可以看到刚才我们节点上其实有三乘，3x2呃，横向三个节点，纵向有三个节点，每个节点两个component是18个自由度，那么polypc，就说我例子上面能不能也用18个自由度。

我例子上面也用18个自由度，那么我从呃假设我只有一个例子，那我从great to的传到粒子上面，就是两个都是18个自由度，如果他们这18个自由度还都是线性相，线性无关的话。

那么我是不是就可以呃做到一个无损的传输，那无损的传输过来，它这个能量的宝石是不是就更好了对吧，这个就是polypc的核心想法，他其实在epic上面更上一层楼，我用更多的mod去呃。

做这个transfer，那这个polyfic，它其实你可以看到，它除了有我们刚才说的这个linear，这个constant和linear这些mod以外，它还有blinear和quatic。

这个实现起来就稍微更复杂一些了，但是它的效果其实也是不错的，这个是一个more more tex的一个模拟呃，你可以看到polypc，它其实就能保持更多的这个mortastic。

然后这边还有一个flip flip，我们会后面再讲，那么flip想法是什么呢，flip他的想法就是我不gather物理量，我们去gather物理量在规的操作上面的增量。

然后我把这个增量gather到particle上面，在gdp的时候，我们只去收集它增量，然后直接把这个增量呃，apply在例子上面，这个增量是啥呢，在你去不可压缩流体的时候。

你增量一般就是这个pressure projection和边界条件对吧，那你做固体模拟的时候，你的增量的就是它裁掉的内里，他的这个内力当然除了内力以外，还有这个边界条件对吧，然后有一些做vfx的同学。

可能会直接用flip来指代一整个solver，那就他可能会觉得flip是一个呃，使用了quin style pressure projection，加上flip。

the vection free ser呃，这个在业界有的时候也会这么说，但是flip它本身只是一个用来遗留的格式，它并不是一个整个的solar，但是有的时候大家会喜欢用flip来说，是整个sol。

这边可能有一点会让人摸不着头脑，所以我额外来讲一讲，这个flip其实也有些历史了，这个flip他是一开始是在jcp，1986年他就提出来了，然后第一次被引入到graphic里面，是2005年。

那么pick flip，他在实现的时候到底有啥区别呢，其实啊区别很小，很pick里其实就是直接从你的呃，下一个时间步的速度场上面gather回来，flip呢他不我刚才说了flip，他不gather呃。

速路上，注意我这边下标是i的话，就是grades，然后下标是p的话就是particle，我去gather什么，我去gather的是新的速度场和旧的速度场的差，我把它gather回来。

然后呢我把它叠加到例子上面，旧的纸上面，那么其实我信息就不完全是从particle degrade，然后归咎particle，对吧，我多了一条particle to。

particle的这个这么个信息路径，这个就能避免很多信息的丢失，如果你纯party to grade，然后graded particle，那你会丢掉很多信息，ok那么flip有什么问题。

flip它就是一般来说有的时候会非常noisy，然后他会看起来噪声非常大，那么pig又过于dissipative，那么有没有什么办法能够取长补短呢，当然有办法啊，这个时候大家就会搞一个叫做flip，0。

99的搞法，什么叫flip 0。99，flip，0。99，就是0。99乘上flip的这个这个操作，加上0。1乘上pc的操作，那你就可以呃，相对来说给flip也引入了一些pick的单品。

pick damping是非常非常大的，你就算只有0。10。01的pick，那你这个单品已经足够了，我们来看看demo，你可以看到，嗯这个ok这个地方有一个flip blending。

这个显示上稍微有点问题，我就如果说我把它调成一，就是纯flip的话，你可以看到它表面是非常noisy，它这个例子就是整个，你可以看到它表面非常的buy对吧，他这个效果不是很好。

这个flip他其实也有很多种实现的方法，很大的一个差别是你做完速度更新以后，你到底是用pick的速度场去移动你的速度，还是用flip的速度场去移动你的速度，这个实现呢是用flip速度上去移动你速度。

它实现出来就会非常非常的easy呃，如果你用pk速度上去移动，速度还会就会稍微好一点，但是这个里面就没有实现另外一种方法，那么我们怎么样让它变得不挪isy的一种搞法，就是你用pg对吧，我把它移动移到零。

你看我用纯pick，那它就完全不同a z，但它也非常的不energetic，他就没有能量，能量很快就耗散掉了，那么我就希望去把它中间来个差值，pick和flip的差值，你可以看到。

这样的话他就没有那么noisy，然后它的也能保持很多能量，这个就其实比较好，一般来说大家呃在用flip的时候，都会blend进去一点，pick这个demo大家可以去玩一玩。

ok那么讲了这么多格式啊，那到底用了用哪一个呢，那么我的建议呢就是如果你刚开始的话，建议你使用epic，为什么呢，因为pc基本上不用pick，它过于呃过于单品了，他就是单品过的。

他这个能量很多模式很快就没了，然后为什么不建议开始用popolypc，或者用flip呢，因为polypc商相对来说实现起来更复杂一些，那么如果要精益求精，你完全可以用polypick。

但是它实现起来稍微稍微复杂一些，然后flip flip的话就实现起来，然后它实现起来会比较noisy，它比较noisy，而且它的实现起来的时候。

你得多开一个grade来去存新的grade和旧的grade，这样你才能算出grade dela对吧，你就要多allocate一个great呃，这个各种feel，所以实现起来稍微复杂一些。

刚才演示的只要三行代码就可以实现了，对于存储非常友好，他不需要去备份你的速度上，同时呢它保持你的angular momentum，看你的角动量是恒的，而且非常stable，然后视觉效果也比较好。

相当于他在各种方法里面站在一个呃，甜点是swiss boss这么一个位置，后面我们会讲m m p m是在干啥哦，有同学有一个问题是particle特别多，内存还是比flip多，他我觉得他意思是。

当你粒子占的内存很多的时候，那网格上的内存是不是就可以忽略不计，那这个完全是讲呃，有的时候粒子上占的内存会很多呃，特别在m p m里面，你的例子会有很多，m p m也是一种混合方法，后面我们会听到。

然后你如果是在模拟流体的时候，你的每一个great cell里面，二维里面你往往会有至少四个例子，三维里面至少有八个例子，所以这样看起来的话，其实例子上用的占用的空间啊，是比规则上面占用空间要多很多的。

所以确实很多时候规则的空间。

大家是忽略不计的，这个是个很好的问题，那么我们就讲到了npm，终于讲到了npm，这个课最终极的目标其实就是讲npm呃，m p m它是一个混合欧拉朗日模拟的scheme。

之前我们讲scheme嘛都是只有a vection对吧，我们呃你要做projection的，你得用一个别的multique server，或者m g p c g to solve，它的压强。

但是m p m它就是一个stand alone的一个schema，他，拉格尔例子和欧拉网格来完成整个模拟的过程，然后如果既然它是一个整个的simulation，scheme的话。

那你particle它存储的东西就不只是速度，它要存好多好多东西，比如说前面梯度，比如说你的现在的体积，比如说呃，你如果是要模拟一些各种别的材料的话，你可能还得存各种别的物理量。

npm在strap里面是非常呃热门的话题，最近几年是非常火，2020年可能就有至少五篇paper是桨叶片，然后最早的时候他是south ky他们去研究出来的，然后在1996年左右。

最早被大概这个时候被提出来，然后第一次引入到graphics里面，然后在graphic里面发扬光大的是，第一次引入的是2013年，然后最早的时候是用来模拟血，最近几年m p m是非常非常热门的话题。

有好多好多paper都是gb m p m，这个polo上面有很多可以研究的，比如说你polo上面各种材料模型，你可以去研究，可以研究血，然后泡沫，然后沙子啊。

各种各种各样grade上面呢你用规律用什么data structure，可以用space grade，可以用open bdp，你可以用multiple grid，多个grid。

然后你当然也可以用太极来表示它的grid，然后particle degree transfer，这中间transfer刚才也讲了很多，比如说pick a pic，然后polypc还有gmp对吧。

还有moving these quares，后面我们要讲m s m p m就是moving this quest，transfer，还有一个叫compatible fat的格式。

它是其实是用来模拟切割的，以及m p m和钢铁的耦合，这个课程里面我们就不讲cpu。

但是m s是一个非常有意思的东西，我们会讲到刚才说到呀，呃最早的时候m p m是按模拟学的，然后这个是迪士尼的这个电影叫冰雪奇缘，然后当时npm还是非常慢的，这个场景模拟了。

据说是模拟了一个星期才模拟出来，当然现在大家有gpu的大家，然后有m s这些各种技巧以后，m p m其实已经可以跑的非常快了，然后再加上太极这个编程语言的加成，所以现在大家可以基本上用88行代码。

写一个能在自己电脑上跑的npm，基本上没有任何难度呃，可以看到它这个整个算法的发展过程，它的复杂程度和性能都是在不停的在被优化的。

mm到底为什么poplar它有很多很多很好的性质，比如说它能自动的把不同的材料给它混合起来，你可以在m p m里面模拟这个立体模拟固体，然后关键的material这个叫颗粒状材料。

比如说沙子或者米饭什么的，都可以用来模拟，然后呢它还可以自动handle collation，自动的字体和字体的碰撞，然后别不同材料之间的碰撞都可以用，然后他的它会有自动的fraction啊。

你一个m p m，你把它拉伸拉的太长了以后，它自动会就会断裂，但是这个fracture并不是数值，并不是一个呃准确的反馈，但是很多时候这个fire也够用了，后面有时间的话，我们可能会提到他有m p m。

也可以和这个连续伤害模型cdm给它整合起来，然后有一套方法叫cd m p，是一个最近的一个比较好的saber paper，后面如果有时间我们也会提到m p m呢，还非常擅长模拟大型变。

这个大型变在graphics里面做，visual effects里面是非常受欢迎的，对吧，你这个一个角色，你总是希望他这个扭来扭去，或者耳朵拖得非常长，或者呃被某一个子弹打了以后发生很大的变形。

就很夸张嘛，大家很喜欢看这样的效果，在动画里面非常有用，然后他既用了呃格子应用的例子，它是一个hybrid欧拉拉格朗日方法，这个欧拉网格啊。

在这边除了其实这边就不是用来做pressure projection，但你也可以做pressure projection，但更多时候大家用欧拉网格是为了去resolve collection。

self collision，然后两个粒子碰撞与碰撞的过程。

是通过欧拉网格来实现，但是大部分物理量还是存在粒子上面，ok那么我们来看一看最早的时候，npm他的这个simulation cycle是好有好多好多部，他有接近啊，不是接近他，就是有十部。

它是首先你有state，然后你做p图g到grade上面，然后grade上面呢你还得去呃asmate下particle volume，这一步我们现在已经不做了，然后你在龟上面还得算这个算它的受力。

然后各种更新，然后帮边界条件，有可能你还要做implies solver，你得做一个影视的时间积分，这个我们下一讲会讲，然后更新到例子上面，然后这样不断的一个流程，可以说还是比较复杂的，你要实现起来。

当时聊实现m p m还是相当复杂的，我记得我当时实现的时候，可能实现了有个有个一个月，那个时候我还是个本科生，然后什么都不懂这个当然。

但是这个算法本身就比较复杂，然后这个是他的早期的时候的实现方法呃，这个大家可以去看一看，还是思路比较清晰的，但是步骤比较多，但是今天我们就不讲这个复杂的方法，我们直接讲简单的讲。

我们在2018年的时候提出来这个mm p movie girl，s p，所以呃他是一个其实是一个比较新的方法，因为它的前前置工作，而是在这个科技树的前一个节点，在16年才解锁。

然后我们在2018年是我本科的时候，在蒋老师实验室，可能是一个暑假做的一个一个项目，还是挺有意思的，然后他有什么好处，他需要用的flops就是浮点数操作，是比呃传统m p m要少了1/2。

相当于他快了两倍，因为mp m它是一个computer b操作，所以我们还是希望去降低它的负电数次数的，如果是mei bb，降低复联数也没有意义，没有意义，但它是个computer b。

然后它比传统mp m容易实现很多很多，非常非常容易实现，然后用太极的话，88行代码就可以实现，然后太极里面有好多demo，大家可以跑一跑，我可以给大家演示几个。

今天呢我们会后面会go through一下，这个m88 是怎么实现的，我这边就给大家呃跑一跑这个m p m88 ，m88 ，其实就是一个最简单的m p m fluid solver，对吧。

它是一个比较简单的，只有一个材料，再看一下它是其实效果还是不错的，你用了apex，就它能比较好的保持面的这些窝的这些效果，然后呢太机里面还有一些别的例子，比如说npm 99，那我就直接把这个拖到这了。

pm 99，99行代码呢，这个里面就有三种材料，有一个这个液体，然后有一个红色的，这个是呃果冻，你可以认为是弹性物体吧，各种都可以，然后嗯白色的是雪，白色的是血，但你如果觉得这个不可交互。

你可以跑m p m21 128，他其实每个代码都差不多，只是加了几条代码，就是更多的材料以及一些用户界面之类的，这个里面你可以用通过拖动鼠标啊，来去操作这些mp例子，这个好像按r k reset。

如果把它不小心全部搞坏了，你可以按一下r reset就可以了，这个是没有重力的，我可以按一下w s d给它加重力，啊反正就是一个几比几几比较好玩的一个demo。

大家可以去跑一跑。

那么接下来我们讲一讲，这个88和代码到底咋实现啊，其实npm啊，和之前讲的apex其实已经非常非常像，m s m p m基本上就是epic，再加上十几行代码去加入npm的一部分。

然后他和之前的party excel啊，这些hybrid alering，a ring的方法其实是非常接近的，你首先要做一个party degree，这个呃把粒子传输到格子上面对吧。

然后再做great operation，然后再做great to particle，其实和之前一模一样，唯一的区别呢就是说你在做p6 g的时候，你可能得算一算每个例子。

它的coaches stress以及它的coaches stress，对于每一个greek note上面的贡献好。

那么我们来看一看具体代码呃，首先我们看这个88行代码的前21行，那其实就是一个很典型的太极代码的前半部分，他前半部分一般是我们去定义一些常数对吧，我们有8192个例子好，当然最重要的请记得要初始化太极。

得指定gpu后端，然后他会在coda或者open gl啊，或者在metal上面运行，然后我们的grade呢，是128x128的，ok然后我们指定一个时间不长dt。

然后每个party呢有一个自己的boom和density，然后质量呢就是volume chandacity对吧，然后我们还有一个bug modules，这边e它就是它的sten，在液体里面就是bug。

modulus，体积模量，ok然后我们来看一看他的subs step kernel，substep knel里面第一部分就是p6 g对吧，你会发现它非常非常像，唯一一个区别就是我多了一个stress。

这个stress是什么呢，stress就是液体被压缩以后，它有一个抵抗压缩的过程，这个它就会在通过一个coaching stress来体现出来，它就会给每一个粒子都会去推开周围的例子。

当然它不能m p m的例子和例子，从来不直接用的，它是通过先推开周围的网格来实现这个操作，那对于这个stress呢，你看我这边除了呃我这边还有刚才还是一样。

epic就是我把每一个例子的transitional velocity，切给它，叠加到grade node上面，然后还有个i fine这这么一个成分，呃，fine里面刚才apex里面的fn就是c。

但是m s m p m里面的f又多，除了c以外还多了一个stress的分量啊，这也要稍微注意一下，它有一个stress分量，那么其实你可以看到m m s m t m是高度，做做最小化改动一个好处是什么。

第一它省计算，因为我可以我这边可以把stress和epic，他这个c给它叠加到一起，统一的用一个fm matrix来来存储，然后这样我在scatter的时候，我就不需要先scattered它的力。

再scatter它的呃，这样我的计算量就少了一倍，然后还更简单对吧，这个地方其实你去看传统m p m，这个地方实现起来还是不容易的，那么接下来啊，great normalization。

这个大家应该都很熟悉了，刚才讲到，我们只需要对质量大于零的地方去做，normalization，这个44行前面其实就是pc epic，polypc前面的normalization，45行以后是啥。

45行以后是边界条件对吧，如果说一个节点它靠近边界太近了，我就得去给把它的速度的某个分量给它设成零，这是最简单的边缘条件，大家可以去理解一下到底在干啥，这边由于时间关系，我就不细讲了。

最后一步是great parle，greg particle和epic其实也非常非常像吧，那其实这边基本上就是完全就是epic，68号之前都完全就是epic，唯一一个区别就是71行，71行。

我们这边更新了一下j j是什么，j是这个粒子的体积，我们怎么更新j呢，一加上dt乘上epic，它的这个c的trace，epic的c的trace是什么，我们来回忆一下a p h。

这个c矩阵到底是什么物理含义，你看到c矩阵还有四项，其中有两项，这个c00 和c11 c0 c一是dation，它的分量对吧，这个dation是c里面，唯一会影响体积变化的两个分量。

因为你如果去share你的材料的话，它体积是不会变化的，当然你只有用这个c00 ，它的你去在横向上去拉他或者c1 ，在纵向上去拉它，它的体积才会变化，然后cd加c一是什么，就是c的trace。

那么c的trace乘上dt，再加上一就是你的呃体积的变化率，那么我们把这个j去乘上这么个体积的变化，就可以得到下下一个实验部的j，那么这边是一个fluid的这个流体的模型。

流体模型里面的力完全就是由j来控制的，那么最后呢就变进到了我们主程序部分，我们得初始化这个例子对吧，然后我们创建一个窗口，在窗口里面不断的去调这个subset，注意你每次做m p m这个cycle之前。

你得去把你的grade上面的速度和质量，都给它设成零，所以啊这就是88行代码哈，其实只有87行代码，但是为了中国人喜欢八嘛，就是我们就叫它88行代码版本就ok了，这88行代码版本还是比较容易实现的。

这个大家每个人应该都可以去尝试实现一下，好那么我们来回顾一下，今天呢我们主要讲了，混合欧拉拉格朗日的这些scheme，那么他的核心思想呢，就是用particle来track你的material去追踪你。

material，你的可能是material上面可能会存位置，可能会存速度，如果你是在做npm的话，你听说上面还会存在formation，用rider呢来去更新你的力的场或者速度的场。

你可能在fluid里面，你会在规则上面做一下quin projection，在你的npm里面，你可以在规则上面去算一下koshi stress，在规则上面的导致了这个内力。

这样我们就有了一个pick这个方法，然后为了减轻pg dissipation，有两套模式对吧，第一种是我们传更多的，传更多的这个information，我们可以加入更多的mod，比如说epic。

epic加入了一个局部的iphine速路上，polypg呢，polpg加入更多的局部的polynomial速度上，另外一种方式呢是用flip flip，不去传输我们的这个所有的物理量。

我们去传输物理量调查对吧，这样就变成了flip这个方法，当然pick和flip可以blender在一起，最后呢我们简单讲了一下，它增加了代表量其实非常非常少。

你是相当于是用particle去存你的deformation information，然后再用epic去传输，你的例子上的各种受力也好，deformation green update也好。

这个理论我们下一讲会提到，大家只要记得m s m p m和epic，真的是非常非常相近了，如果大家想学习更多npm的内容呢，大家可以去看张老师和joe特曼，还有angel sally这些。

其实呃npm最早的五个researcher，他们在sf 2016年上做的这个课程，还是质量非常高的，他们有个cosnow是cosnow的。

堪称呃最我看过的最well written the course note，它整个公式啊都非常的清晰，大家可以去看一看，其实里面也介绍了各种连续介质力学啊，可以说图文并茂。

学习起来的事还是非常非常不错了。

我们最近2019年，我们在surah上面还有一个呃，如果大家对高性能感兴趣，可以去看一看这个高性能的cosnote，然后如果说大家想学更多npm的paper的话，可以看蒋老师的paper list。

这个是一个npm最近的paper，所有的配合者集合在蒋老师的主页上，ok那么今天就到这了，我看看有没有什么问题需要解答的，还是有不少问题的，啊我有同学说p to g的一部分，这个加等于。

因为你不同的例子，可能会同时去accumulate到同一个good node对吧，这边就有个data race，那么这个加等于在太极里面是自动的，omic，这个我们之前其实也提到过呃。

有的版本里面可能会写tian atomic ad，把它显示的写成atomic，但我们这边是用的是自动atomic，在太极里面所有的加等于都是自动的tommy，如果你不要加，不用自动的tommy怎么办。

你就写成规的v，这个等于规律v加上它其实就你不用加，等于你用a等于a加b把它显示写出来，那就是butommy怎么做到不可压缩的啊。

这个是个好问题，这个例子真的不可压缩吗，这个m p m128 88，这个里面它其实是弱可压缩的一个weekly compressible，它不完全是不可压缩的，这边其实是用他的。

你可以看到我刚才提到bug mojs对吧，其实这个是一个可压缩流体啊，但是它看起来不太可压缩，就是这个其实它是可压缩的，那它怎么保持体积的呢。

是因为我们通过维持了维护了这个j，然后这个j你可以看到它，j这边会算出一个stress，他这个stress是jp减去一，然后乘上了一坨系数对吧，呃当你的j小于一的时候。

它就会有一个推开周围格点上速度的这个题，至大于一的时候，它就会有嗯拉近周围割裂上速度的这么一个力，所以它其实是通过j来保持自己的体积。