这样也可以 ， 看来我还是得看看 cpp11
# case 3 test okay
for i, j, k in ti.ndrange(grid_n, grid_n, grid_n):
# we cannot visit prefix_sum[i,j] in this loop
pre = ti.atomic_add(prefix_sum[i,j], grain_count[i, j, k])
linear_idx = i * grid_n * grid_n + j * grid_n + k
list_head[linear_idx] = pre
list_cur[linear_idx] = list_head[linear_idx]
# only pre pointer is useable
list_tail[linear_idx] = pre + grain_count[i, j, k]
- 其实文中提到的并行计算方式完全可以用 ti.atomic_add 去掉所有串行区域，进一步并行化，你觉得应该咋搞？
对应这个问题 提了一个PR给你
请问在你的代码中折半邻居搜索后计算效率会有所提升嘛，我试过折半搜索，同时去掉了相邻网格时id1要小于id2的条件，照理来说减少搜索次数速度应该提升，但发现效率反而下降了
会提升的，i<j的条件需要选择性保留
《偷功》
实现了一个简单旋转小球的效果
效果：
PBF 版本： 增加了不可压缩性，粘度等， 这个实现参考了 Ten Minute Physics 中的实现
代码：GitHub - mrzhuzhe/taichi_dem: A minimal DEM simulation demo written in Taichi.
python3 dem-3d.py # 黄色 dem 旋转小球
python3 pbf.py # 蓝色 pbf 旋转小球
python3 dem-3d-snode.py # 黄色 dem 旋转小球，但是临域搜索的链表结构是用dynamic snode实现
已实现：
- 事实上我们可以将邻域搜索的范围从 3 x 3 进一步缩小一半左右，你能想想怎么构建这个搜索算法吗？
- 其实 Taichi 对于 GPU 上的动态数据结构有很好的支持，如 taichi_elements 中就采用了 dynamic SNode 来维护每个网格 cell 中的粒子下标数组。你能尝试用 Taichi 的dynamic SNode 来实现本文的功能吗？
- 其实 Position-based fluid (PBF) 也可以用这种方法进行加速，你可以试一试吗？
- 其实文中提到的并行计算方式完全可以用 ti.atomic_add 去掉所有串行区域，进一步并行化，你觉得应该咋搞？
后续问题：
if (pos[i][0] >= maxX):
pos[i][0] = maxX - eps * ti.random();
岩土专业，代码小白。在实现过程中遇到一点问题，想要请教社区的大家
代码在这：GitHub - chenqing253/taichiDEM
尝试实现加入tangential force and rotation, 遇到了以下几点问题：
t_i is Tangential unit vector. 上述就是求两个颗粒在接触切平面上的相对滑移来求摩擦力
可以看出切向力可能需要在前一步的接触上的切向力进行更新，有什么方法或者数据结构可以记录 前一步的接触吗？
我用红色和蓝色标注了不同旋转方向的颗粒，旋转了但是低于平均值的为白色。而且颗粒模拟过程中比较Q弹，我也想不出哪里出错了，是因为一定要在前一步的切向力上跟新吗？
赵仕威老师有篇文章里有写A thread-block-wise computational framework for large-scale hierarchical continuum-discrete modeling of granular media，采用多个列表存储数据
我也来分享一个小的demo~，我在代码中加入了切向力，即摩擦力，实现方法比较简单粗暴。为了实现静摩擦，添加了两个integer类型的ti.field容器用于‘记忆’上一个时间步和现在的接触，接触用两个接触颗粒的id的数对表示，还有两个float类型的ti.field 容器用于储存相应的切向力。Field容器的形状为nx10, n 为颗粒的个数，10是考虑到这个颗粒周围的接触不会超过10个。求算接触 (i,j) 的力时，会遍历第i行的容器。（本来想直接在class里加一个list的…)
第一次尝试用taichi语言，实现过程中踩了不少坑，例如调参数总是把试样调爆炸 orz 。希望随着不断学习taichi的过程中能实现得更漂亮简洁一些，像用 n \times n 的稀疏数组来实现对接触的记录。
代码仓库：A small demo implement shear force based on Taichi DEM
代码效果：
加入切向力后可以‘支棱’起来：
没加入之前总是会变得比较平：
这里没有放gif，因为我调小了重力整个模拟过程很慢，一点也不炫酷orz，否则加入或不加入切向力的效果是一样的，我尝试调过stiffness，但是会‘爆炸’，以后还得再调调参数。
参考文献：
Bell, N., et al. (2005). Particle-based simulation of granular materials. Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation. Los Angeles, California, Association for Computing Machinery: 77–86.
Lee, J. and H. J. Herrmann (1993). “Angle of repose and angle of marginal stability: molecular dynamics of granular particles.” Journal of Physics A: Mathematical and General 26(2): 373-383.
从工程角度使用太极语言完整实现离散单元法。
使用BIMBase ™二次开发可视化。BIMBase ™是北京构力科技有限公司研发的用于BIM的图形平台。BIMBase开发平台-BIMBase Platform
使用Houdini ™可视化。Houdini ™是一款特效软件，工作流程灵活，对用户友好。Houdini | 3D Procedural Software for Film, TV & Gamedev | SideFX
因Taichi论坛上传图片大小限制，上述Demo GIF效果欠佳。请前往GitHub代码仓库查阅高清GIF (地址下述)。
Denver Pilphis (Di Peng) - 离散单元法理论与实现
MuGdxy (Xinyu Lu) - 性能优化
提交分支 – 按照比赛规则要求，去除对其他库的import引用，可能导致运行效率变低：
开发分支 – 持续迭代和优化：
本例提供一个完整的离散单元法 (DEM) 仿真的实现，考虑了复杂的DEM力学机制，仿真结果达到工程定量精度。
本例使用Taichi语言，加上适当的数据结构和算法，以保障计算效率。
相比Taichi DEM原始版本，本例增加了如下功能：
二维离散单元法→三维离散单元法；
完整考虑并实现颗粒方位和旋转，预留了建模和仿真非球形颗粒的空间；
实现了墙 (离散元中的边界) 单元，并实现了颗粒-墙接触解算；
实现了复杂离散元接触模型，包括一个胶结模型 (爱丁堡胶结颗粒模型，EBPM) 和一个颗粒接触模型 (Hertz-Mindlin接触模型)；
因胶结模型已实现，可以使用胶结团块模拟非球形颗粒；
因胶结模型已实现，可以仿真颗粒破碎过程；
将材料特性与颗粒/墙体关联，以降低空间成本；
将接触面交互属性与接触关联，以降低空间成本；
基于Morton编码实现空间哈希表，进行邻接搜索 (宽碰撞检测)；
分别对低工作负载和高工作负载采用位表和并行扫描算法，存储邻接对以减少内核内部的发散，从而提高并行计算效率；
接触存储在每个颗粒所关联的动态列表中，以降低空间成本，列表在每个时间步都会维护（处理接触的增减）。
原理将另行上传附件说明。
本例展示了开仑台球开局第一球。白球正对其他球运动并发生碰撞，然后球散开。虽然该过程中有能量损失，但是由于没有抗转动力学响应耗散转动动能，所有的球在进入纯滚动状态后将会一直滚动下去。本例可用于验证Hertz-Mindlin模型。
本例展示了一个立方体胶结团块撞击平整表面的过程。
撞击过程中，团块内部的胶结受力将发生破坏，然后团块碎成碎片，飞向周围空间。
本例可用于验证EBPM模型。
本例与上例相似，唯独团块所含颗粒数不同。本例可用于大规模体系仿真的性能测试。
本例展示了一个斯坦福兔子形状的胶结团块自由落体撞击水平表面，然后破碎的过程。本例可用于大规模体系仿真的性能测试。
感谢谢菲尔德大学化学与生物工程学院Xizhong Chen博士对本研究给予的帮助和支持。
附件补充原理说明：
HPC_Part.pdf (410.2 KB)
DEM_Part.pdf (307.9 KB)
此外，因为发现不同学科对于同一概念使用的术语不同，为避免歧义，我们约定统一使用以下术语，这里进行补充解释。
DEM一个时间步包含5个流程：
邻接搜索 - 接触判断 - 接触解析 - 接触评估 - 时间积分
Neighboring search - Contact detection - Contact resolution - Contact evaluation - Time integration
邻接搜索 Neighboring search == 碰撞检测宽检测 Broad phase collision detection
此过程解决两个颗粒 (==粒子，英文均为particle) 是否相邻，因为只有相邻的颗粒才有可能接触。虽然该过程对于DEM算法本身并非必要，但是此过程可以剪枝很大一部分冗余计算，大大提升计算效率。
接触判断 Contact detection == 碰撞检测窄检测 Narrow phase collision detection
此过程解决两个颗粒是否真正有接触 (== 碰撞, contact == collision)。接触的意思是，二者之间存在一个相互作用 (interaction)，可在颗粒间传导力 (矩)。
接触解析 Contact resolution ∈ 碰撞后的物理计算
此过程解决两个颗粒的接触的几何特性，以作为接触评估时输入接触模型力-位移表达式的位移数据。
接触评估 Contact evaluation ∈ 碰撞后的物理计算
此过程调用预设的接触模型 (contact model)，接触模型中含有力-位移表达式来表达接触的力学响应。通过输入位移数据，计算后得到接触传导的力 (矩)，然后作用于接触两端的颗粒上。
时间积分 Time integration
此过程汇总颗粒所受的所有力 (矩)，利用牛顿第二定律和时间积分得到颗粒的运动学数据，然后更新体系。
邻接对 Neighboring pair
为两个经邻接搜索确认相邻的颗粒的编号 (ID) 所组成的数对 {ID_i, ID_j}。
邻接对列表 Neighboring pair list
存储所有邻接对的列表。
接触列表 Contact list
存储所有接触的列表。接触也是以 {ID_i, ID_j} 的模式进行存储，与邻接对类似。
哇，这个投稿好棒呀！！
来报个到
，和彭老师合作特别顺利！短短半个月时间就做出了这样一个有趣的demo。
是22天 (纠正)
确实，按照你补充的网格信息，原文档是说的通的。
更改和补充Demo：
本例展示了一个斯坦福兔子形状的胶结团块自由落体撞击水平表面，然后破碎的过程。本例可用于大规模体系仿真的性能测试。
本例展示了一个斯坦福兔子形状的胶结团块自由落体撞击水平表面的过程。由于胶结强度极高，兔子不会破碎；相反，由于胶结弹性模量相对较低，兔子会表现出相当软的力学响应。本例可作为上例很好的对比。
因Taichi论坛上传图片大小限制，上述Demo GIF效果欠佳。请前往GitHub代码仓库查阅高清GIF (地址详见原帖)。
最后几天截止时间我也来加一点：
附上我的源码以及动画：
这是不加rolling的自然休止角
由于文件大小限制，加了rolling和twsiting之后的自然休止角动图请移步
并且通过计算几个案例的计算，验证了代码的可靠性。
两位如果使用过程中有遇到什么磕磕绊绊的地方，欢迎到 GitHub 开 issue，我们尽力让这个过程更加顺滑 
