粒子与网格归属问题

[夏雨天]

对粒子并行，每个cell生成concurrent vector，用粒子坐标计算出所在网格，pushback即可。

[Number44]

通过findCell()函数获取粒子所在网格单元编号，用Foam::labelListList&以cellLabel-particalLabel的形式存储，应该一次partical循环就可以的。
另外，感觉这事和sdfibm这个第三方库有点类似，这库也要判断粒子是否包含网格，和这事刚好相反，或许可以参考一下。

[北斗]

这题我会！
可以建立两个数组（实际上是两个哈希表）。
一个是以粒子编号为下标，存储的是网格编号。一个是以网格编号为下标，存储的是粒子编号（因为一个网格可以有很多例子，所以这是个二维不等长数组，第二维可以用std::vector）。

每次粒子移动以后，更新这两个数组。

假如知道粒子编号想找所在网格，那个就第一个数组；假如知道网格编号想知道它里面有哪些网格，就用第二个数组。

[同学博]

@马乔

[北斗]

另外，假如你的问题只是邻域搜索的话，即
知道一个粒子编号，想求它所在的网格的其他粒子的编号，
那么还可以用最小堆实现。

PS：最小堆是一个自动按照从小到大排列的数组。（可以在每个时间步开始的时候排序）

我们建立这样一个数组：
这个数组的下标为粒子编号，而存储的值是网格编号。这个数组总是按照网格编号进行排序。
（我们称之为数组a吧。网格编号记作CID，粒子编号记作PID）

由于这个数组是从小到大排序好的，所以它应该大概是这样的

24 35 2 31 42 57
0 0 0 1 1 1
上面的代表粒子编号即下标，下面的是网格编号即值。
总之粒子编号是乱的但网格编号是排好的

假如我知道某个粒子的编号是31，那么直接通过索引粒子编号就知道它在1号网格，即

a[31]=1

那么搜索所在网格的其他粒子的时候呢，只需要同时向前和向后搜索就行了。直到搜索得到的
网格编号和当前粒子的网格编号不等为止，即CID！=1即可。

这就避免了搜索所有的粒子。

[北斗]

这个说的有点问题
排序之后不是下标乱了而是值乱了。
所以其实应该不能用数组。
大概应该用哈希表，键为PID 值为CID。
按照值排序，按照键直接索引得到值。

[李东岳]

感谢各位大佬，太神奇了。目前我在想或许DynamicList或者hash table是个好出发点

[李东岳]

我现在用的最笨的方法弄得：

1. 网格遍历
2. 粒子遍历
3. 如果粒子恰好在网格内，就把label存储在当前网格的DynamicList

        forAll(U_, cell)                                         
        {
            DynamicList<label> pL(0);
                                                                 
            scalar pN = 0;
            forAllConstIter(typename MomentumCloud<CloudType>, *this, iter)
            {                                                    
                const parcelType& p = iter();                    
                if (p.cell() == cell)                            
                {                                                
                    pN += 1;                                     
                    pL.append(p.origId());
                }                                                
            }
            
            if (pN != 0)                                         
            {
                Info<< "Cell[" << cell << "] has " << pN << " particles, "
                    << "particle label is " << pL << nl;         
            }
        }   

输出结果还可以：

Cell[2295] has 2 particles, particle label is 2(20 21)
Cell[2297] has 2 particles, particle label is 2(4 7)
Cell[2298] has 1 particles, particle label is 1(3)
Cell[2340] has 3 particles, particle label is 3(17 18 19)
Cell[2341] has 5 particles, particle label is 5(11 13 14 15 16)
Cell[2342] has 4 particles, particle label is 4(5 6 8 9)
Cell[2343] has 1 particles, particle label is 1(1)
Cell[2386] has 1 particles, particle label is 1(12)
Cell[2387] has 1 particles, particle label is 1(10)
Cell[2388] has 2 particles, particle label is 2(0 2)

但是这里面遍历网格+遍历粒子。肯定是要慢。

[马乔]

@李东岳 这么写呢？

            List<DynamicList<label>> pL(U_.size());
            forAllConstIter(typename MomentumCloud<CloudType>, *this, iter)
            {                                                    
                     const parcelType& p = iter();                                              
                    pL[p.cell()].append(p.origId());                                           
            }            

或者

std::Multimap<label,label> Lp;
forAllConstIter(typename MomentumCloud<CloudType>, *this, iter)
{
     const parcelType& p = iter();  
     Lp.insert(std::pair<label,label>(p.cell(), p.origId()));
}

都是我云的

[李东岳]

@马乔

            List<DynamicList<label>> pL(U_.size());
            forAllConstIter(typename MomentumCloud<CloudType>, *this, iter)
            {                                                    
                     const parcelType& p = iter();                                              
                    pL[p.cell()].append(p.origId());                                           
            }   

这种方法不需要遍历两次。类似一种DynamicListField。不过有很多元素是空的。不过比我的方法要快。挺好。多谢大佬

[马乔]

@同学博 你这么操作map是对的，我开想的是查找key

[同学博]

@马乔

[TM]

这段用来统计网格单元内粒子个数的代码，编译的路径是什么，求东岳老师指导