请选择 进入手机版 | 继续访问电脑版

开源计算机图形学社区(Open Source Computer Graphics Community) |OpenGPU Forum (2007-2013)| OpenGPU Project

 找回密码
 注册
搜索
查看: 1440|回复: 15

写了一个基于fem的柔体模拟 [复制链接]

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-5 11:53:56 |显示全部楼层
终于抽出点时间写了一个简单的柔体模拟器。算是自己几年PHD生涯的一点成果。 分享给大家。希望能有所帮助。一些特性:
1,简单,self-contained. 不用连接什么外接库。
2, Invertibale FEM. 跟其他库只支持C3D4不同,只要实现相应的shape func 和 shape func derivative就可以支持任意solid Finite element. 暂时只写C3D4, 其他的以后慢慢再加。
3, 碰撞检测。 DCD和CCD都支持。
4, 碰撞反应。contact 和 friction 都支持。 实现的是 implicit contact handling of deformable objects. 这个找了好久,都没有找到,就自己写了。

一些代码直接copy的 vega fem, Febio, el topo。这些都是很优秀的开源库。但是找不到 一个完整的能直接用来模拟的库。vega fem 没有碰撞检测和反应。 febio 不支持ccd也太复杂。 el-topo没有simulation。
还有一个 库叫cubica, 很不错。但是碰撞似乎是用的penalty.

还没来得及写comment. 大家先将就看看。以后一些feature 和comment 再加上去。
编译很简单,用qt creator 打开就好。不需要连接什么其他库。一些第三方库如eigen 和 libigl都是header only, 也包含进去了。运行后按空格就可以模拟一个torus 掉在地上。
如果大家觉得有什么feature是特别想加进去的。请告诉我,有时间和能力我会来实现下。

https://github.com/nydragon/zxDeform

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-5 11:58:15 |显示全部楼层
我是在ubuntu上面开发的。也不知道windows下会不会有问题。但是文件组织很简单。所有文件都在src里,根目录里是用qt写的gui, eigen 这些在third-party 里。 直接添加再设定下路径应该就好。

使用道具 举报

Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20

注册时间
2010-3-27
积分
6027
发表于 2017-5-8 09:50:44 |显示全部楼层
欢迎大神参与,能将理论变成产品(代码)。我碰到的大多数计算机“大神”都只喜欢谈理论,见不到它们的代码。
咨询个问题:vega中关于Model Reduction的代码大致在哪个文件?自己找太花时间

使用道具 举报

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-8 15:04:22 |显示全部楼层
ljb 发表于 2017-5-8 09:50
欢迎大神参与,能将理论变成产品(代码)。我碰到的大多数计算机“大神”都只喜欢谈理论,见不到它们的代码 ...

reducedElasticForceModel, reducedForceModel, reducedStvk

使用道具 举报

Rank: 20Rank: 20Rank: 20Rank: 20Rank: 20

注册时间
2010-3-27
积分
6027
发表于 2017-5-10 16:34:53 |显示全部楼层
谢谢,看下再来咨询。

使用道具 举报

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-17 23:31:09 |显示全部楼层
update:
1, 加入了对C3D10的支持,并实现C3D4 到C3D10的转换。
2,实现了基于augmented lagrangian method 的 face to face contact。 这其实是基于febio的一个简化。适合对精度要求较高的continuum mechanics,暂时只支持静态模拟。其实LCP方法不适合用于高精度计算,因为模拟很难收敛,只是速度较快,所以适用于visual effects 或 game。

所有代码都还没有加comments,大家有什么疑问可以直接问我。等大部分feature实现后我会加上comments的。

使用道具 举报

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-20 22:54:36 |显示全部楼层
update:
1, 加入了对六面体(C3D8)支持。
2, 实现基于augmented lagrangian method 对ccd和dcd的contact solve. 这个跟之前的不同,精度要求低,模型线性化简化成了quadratic problem。相比lcp更简单,因为实际上的实现就是一个penalty based method。适合于多种力的模型。
3,加入了reduced model。 这是一个基于cubature的模型。 contact handle 是基于augmented lagrangian method。

使用道具 举报

Rank: 9Rank: 9Rank: 9

注册时间
2012-11-4
积分
352
发表于 2017-5-28 15:09:07 |显示全部楼层
本帖最后由 不要说 于 2017-5-28 15:14 编辑

支持下LZ,我自己也在写基于FEM的固体形变模拟,代码比较挫就不公布github地址了。

就问一下,LZ有研究过unstructured multigrid preconditioner吗?这玩意用来加速椭圆(elliptic)型偏微分方程的求解理论结果似乎挺不错的,可以把矩阵条件数压到对数级别(见[1]),但无奈做动态分析的话其实是要解双曲(hyperbolic)型方程,我目前还没看到有证明到能够压条件数的阶,不过似乎DMM他们试了下似乎对大规模矩阵实际效果还行(见[2]),就是不知道是不是真的这样。
特别提一下,由于自己实现的网格生成器生成的是非结构化(unstructured)网格,所以就不考虑经典的基于结构化网格的multigrid方法了。

[1]Adams, Mark. "Evaluation of three unstructured multigrid methods on 3D finite element problems in solid mechanics." International Journal for Numerical Methods in Engineering 55.5 (2002): 519-534.
[2]Parker, Eric G., and James F. O'Brien. "Real-time deformation and fracture in a game environment." Proceedings of the 2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation. ACM, 2009.

使用道具 举报

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-28 16:44:36 |显示全部楼层
本帖最后由 peterntu 于 2017-5-28 16:47 编辑
不要说 发表于 2017-5-28 15:09
支持下LZ,我自己也在写基于FEM的固体形变模拟,代码比较挫就不公布github地址了。

就问一下,LZ有研究过u ...

这个有过一点了解,但是没有深入进去,因为个人看法是multi-grid在structured grid的表现会比较好,而且这样也可以更加容易并行化( multi-grid 其实严格讲不适合并行化,因为收敛要好的话要用GS iteration。但也不是没有办法。但用在unstructured grid上就更难了)。形变模拟其实是可以看成每一步都在求解一个Ax=b的线性方程组。这个一般是用conjugate gradient来求。multi-grid 也可以看成另一种方法。然而在computer graphics 领域,其实这个求解已经不是瓶颈了。瓶颈在碰撞检测上。对于复杂些的场景,碰撞检测占时超过80%。 对于精度要求很高的continuum mechanics, 其实还是在用direct solver。 因为在continuum mechanics上经常出现材料系数差距很大个柔体。比如人体的肌肉(0.15 Mpa)和骨头(7300 Mpa)。这个如果用迭代算法的话收敛会很慢,甚至会使得整个模拟都不收敛。个人觉得Multi-grid在这种情况可能比cg好点。因为可以用更coarse得网格将high frequency的量去掉。而cg每一步的A都是对相同大小的网格。

个人觉得柔体模拟其实最大的难点就是contact或者叫碰撞。因为这里有很大的非线性成分,线性化很困难。在高精度计算中甚至认为一般的c0连续的三角网格是很不精确的。现在新提出的isogeometry analysis似乎是个有潜力的方法。如果能考虑multi-grid方法在这种contact 方向而不是一般的 non-contact deformation simulation收敛速度方向会很有应用。

使用道具 举报

Rank: 9Rank: 9Rank: 9

注册时间
2012-11-4
积分
352
发表于 2017-5-28 16:59:10 |显示全部楼层
本帖最后由 不要说 于 2017-5-28 17:02 编辑
peterntu 发表于 2017-5-28 16:44
这个有过一点了解,但是没有深入进去,因为个人看法是multi-grid在structured grid的表现会比较好,而且这 ...

先感谢LZ的解答吧,不过LZ回答的似乎和我要问的东西不一致。。。。。。
multi-grid其实可以用来当conjugate gradient的preconditioner,因为conjugate gradient收敛性基本上是用矩阵条件数刻画的(当然,准确来说和矩阵特征值的分布也有关)。用multi-grid来做preconditioner就是想去压低条件数的。不过,还是如前面所说,这个条件数的改进似乎对于双曲型偏微分方程,至少理论上本人没发现什么阶的改进。你这边有实际实验过吗?
另外,其实我也列举了非结构化multi-grid的文献,要求结构化的话感觉还是太受限了。
不过就你个求解不是瓶颈,我也不大认同,假如稍微夸张一点到上万自由度的话,现在的性能也没到30ms吧。
Edit:
isogeometry analysis我个人认为也不错,计算力学的一些大佬也在推,而且似乎不用那么麻烦专门去生成网格了。

使用道具 举报

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-28 18:07:36 |显示全部楼层
不要说 发表于 2017-5-28 16:59
先感谢LZ的解答吧,不过LZ回答的似乎和我要问的东西不一致。。。。。。
multi-grid其实可以用来当conjugat ...

multi-grid 我确实只能说了解一点。只能大概说说,不对的话请见谅。

multi-grid可以看做一个preconditioner也可以看成一个solver。因为他这两个是couple在一起的。如果是看做preconditioner的话跟 cg 的其他preconditioner的比如jacobi preconditioner是一个作用。也就是你说的降低condition number。但是原理是不一样的。要想降低condition number首先就要知道condition number为什么会很大,在哪里会很大。我个人的经验是大的condition number来自四个地方:1, 材料系数差距很大的柔体。很stiff的部分相对很柔软的部分会造成很大的condition number。2, 质量很差的网格,比如一个很扁的四面体就足以使整个系统的condition number变得很大。3,同样一个物体,使用更finer 的mesh的condition number 比coarse mesh 大。4, 在有contact的时候,通常需要用penalty的方法来约束contact。如果这个penalty很大,同样会造成condition number变大。

multi-grid方法其实是集中在fine mesh比coarse mesh的condition number大这一点上。通过在不同的level上求解可以进行在这个level上更有效的处理。这样2,3的问题就在不同程度上得到处理。但是1,4号问题不是很清楚怎么解决,所以大家通常是先用比较小的penalty然后在逐渐增大。其实05年有一篇论文已经采用了非结构化的multi-grid方法A Multigrid Framework for Real-Time Simulation of Deformable Volumes (vriphy)。作者claim还不错。

我不是数学系出身,所以对诸如双曲型偏微分方程之类的就不是很懂了。以上只能说是根据自己的一些理解提出的看法。肯定数学上是不准确的。

至于求解不是瓶颈,我是指的在computer graphics主要应用在vfx和game上而言,这一类应用的碰撞复杂度一般而言是远大于物理方程的复杂度的。尤其是对于衣服布料模拟。大部分情况下的方程condition number不会很高,用cg就好,而且也容易使用gpu加速。上万个自由度其实是很小的模型,3k的点的物体的求解现在很容易达到30ms了。

在continuum mechanics领域,求解的复杂度就远大于求解了。因为这些工程问题对精度要求很高。网格密度非常大。contact surface往往都是预先指定好的。碰撞的检测不是很大的问题,难点集中在碰撞反应的收敛上。这个领域应该使用multi-grid的机会大些。但我个人觉得因为精度要求很高,所以收敛性比速度要重要很多。如果把方向集中在这个方向会更有应用。

isogeometry最开始的motivation确实是不用去麻烦的生产网格。但是我个人认为更大的潜力是这个几何体直接就是高阶连续的。不像三角网格是c0连续,contact的收敛会好很多。复杂边界条件可以精确表示,这个太重要了,因为很多工程问题的精确度都是边界条件决定的。

使用道具 举报

Rank: 9Rank: 9Rank: 9

注册时间
2012-11-4
积分
352
发表于 2017-5-28 18:42:00 |显示全部楼层
本帖最后由 不要说 于 2017-5-28 19:19 编辑
peterntu 发表于 2017-5-28 18:07
multi-grid 我确实只能说了解一点。只能大概说说,不对的话请见谅。

multi-grid可以看做一个preconditio ...

蛤蛤,辛苦LZ说得这么详细。但其实,你应该也可以看出来,我自己也是懂点mulitgrid的了,LZ说的为什么条件数会比较糟糕,表象上的原因我这边看到的也差不多,其实如果材料本身刚度高的话,条件数也会变差。数学的东西我也不再提了,其实为什么对于椭圆方程multigrid表现这么好,而对双曲型方程却很困难,这个可以看出来和方程的性质是有很大关系的。而且对于条件数和网格密度的关系,其实有很经典的论断,LZ也可以看一下。

Edit: 试了下,其实单纯求解的话,上万自由度30ms还是可以的,只是整个流程下来有些够呛。

其实continuum mechanics这个圈子划得也太大,毕竟我们的讨论还要把流体去掉,是吧:)。其实就单说固体力学,超弹性,粘弹性,弹塑性都有各自的东西了,其实像incompressibility和nearly incompressibility做起来就已经挺麻烦的了。

关于isogeometry analysis,其实不用麻烦生成网格还是很关键的。你写一个功能比较完整的3D网格生成器就会知道这玩意有多麻烦,尤其各种dirty的输入怎么处理真的会很恼人。

最后,还是说一下,LZ也不用答得这么辛苦吧,针对问题说说其实就行了。

使用道具 举报

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-5-28 19:48:19 |显示全部楼层
不要说 发表于 2017-5-28 18:42
蛤蛤,辛苦LZ说得这么详细。但其实,你应该也可以看出来,我自己也是懂点mulitgrid的了,LZ说的为什么条件 ...

好吧, 我确实对multi-grid在condition number上的作用不是很懂。以上的说法是只是基于我去年实现一个基于fem的人体脚的模拟总结出的。至于大牛们的经典论文了解不多。 unstructured multi-grid也就了解上面提到的那一篇论文。
至于现在的速度,不知道你说的整个流程是否包括碰撞。如果不包括的话
Projective dynamics: fusing constraint projections for fast simulation claim 是30k的自由度可以达到实时
Quasi-Newton Methods for Real-Time Simulation of Hyperelastic Materials 是最新的论文,还没有看,具体能多快还不清楚。
GPU 加速的有multi-grid也有cg-based。大概100k个点可以达到实时。

使用道具 举报

Rank: 9Rank: 9Rank: 9

注册时间
2012-11-4
积分
352
发表于 2017-5-28 20:08:25 |显示全部楼层
本帖最后由 不要说 于 2017-5-28 20:17 编辑
peterntu 发表于 2017-5-28 19:48
好吧, 我确实对multi-grid在condition number上的作用不是很懂。以上的说法是只是基于我去年实现一个基 ...

我惭愧地表示,3187个nodes,9213dofs,我跑一个完整的步长(只含一次牛顿迭代)需要50-60ms(当然不包括碰撞)。求解器全靠自己撸,看来还是得弄个BLAS来用用了。
关于multi-grid用作preconditioner的话,前两年有篇关于preconditioning的综述讲了不少,不过印象中都是关于椭圆方程或抛物方程的,双曲方程的没有 (Wathen, Andy J. ”Preconditioning." Acta Numerica 24 (2015): 329-376.)
你给的第一篇文章看了一眼,似乎不完全是FEM是吧。

使用道具 举报

Rank: 12Rank: 12Rank: 12

注册时间
2016-3-13
积分
865
发表于 2017-6-2 10:51:14 |显示全部楼层
本帖最后由 peterntu 于 2017-6-2 10:59 编辑
不要说 发表于 2017-5-28 20:08
我惭愧地表示,3187个nodes,9213dofs,我跑一个完整的步长(只含一次牛顿迭代)需要50-60ms(当然不包括 ...

你如果是CPU实现的话。求force和stiffness matrix的时候可以用openmp做并行化。解线性方程组可以先用cg。如果收敛太慢迭代次数太多的话就用direct solver。推荐使用pardiso。这个是mkl里面自带的。现在有免费版本可以用。eigen library实现了接口,很好用。

projective dynamics确实不完全是fem,但他的一种特殊形式是完全的fem。当然也只是对于一种很特殊的材料模型。所以第二篇似乎就好很多。但是还没有看。
或者你可以试一下vega fem 这个库。是Jernej Barbič写的库。里面用了不少优化。你可以运行一下看看速度怎么样。作为一个benchmark。毕竟速度跟很多因素都相关。不能一概而论。

使用道具 举报

Rank: 9Rank: 9Rank: 9

注册时间
2012-11-4
积分
352
发表于 2017-6-2 22:26:23 |显示全部楼层
本帖最后由 不要说 于 2017-6-2 22:29 编辑
peterntu 发表于 2017-6-2 10:51
你如果是CPU实现的话。求force和stiffness matrix的时候可以用openmp做并行化。解线性方程组可以先用cg。 ...

老实说,确实现在的solver没作多少优化,优化主要也在单核上,即没上多线程也没上SIMD。主要是这是我业余项目,只能时断时续地做一下。:(  solver这块我自己还是倾向于自己先踩一遍坑,毕竟从我自己的视角这块性能还是很重要的,调库什么的后面再谈。
Vega的话本人其实看过,感觉他们solver这一块似乎也没怎么发力,不过不知道现在的版本怎样。倒是DMM讲了一些优化的trick,看上去还是挺有借鉴意义的(虽然DMM的源码就甭想了)。



使用道具 举报

最近看过此主题的会员

您需要登录后才可以回帖 登录 | 注册

‹‹
我的工具栏

关于我们|手机版|Archiver|开源计算机图形学社区(Open Source Computer Graphics Community) | OpenGPU Project | OpenGPU Forum (2007-2013)

GMT+8, 2017-9-20 20:45 , Processed in 0.084797 second(s), 11 queries .

Powered by Discuz! X2

© 2001-2011 Comsenz Inc.

回顶部