粒子与网格归属问题

李东岳

假定粒子的编号是1、2、3、4、5、6、7。每个粒子归属某一个网格单元。比如

粒子编号	网格编号
1	23
2	23
3	26
4	25
5	26
6	23
7	20

也就是说编号23的网格里面，有3个粒子，分别是1、2、6。我在想怎么写一个简单的代码，把每个网格单元的粒子编号提取出来呢？

网格23的粒子编号是1、2、6
网格26的粒子编号是3、5
网格25的粒子编号是4
网格20的粒子编号是7

一种方式是对网格进行遍历，比如遍历到23网格的时候，遍历粒子，看哪些粒子在这个23网格里面。但太慢了。

遍历粒子的话，是不是更好点，就不需要遍历网格了？但这个算法没想出来是不是岁数太大咯

anubis

@李东岳在粒子与网格归属问题中说：

遍历粒子的话，是不是更好点，就不需要遍历网格了？但这个算法没想出来是不是岁数太大咯

假设有m个网格，那就建m个空数组，知道的是粒子对应的网格编号嘛？那就遍历粒子的时候把粒子的序号添加到对应的数组里就可以了吧。
不过这个方法要建m个数组，感觉也不是太好，肯定有更好的办法，不过编程能力有限只能想到这个。

fubianhanshu

@李东岳这个实现了就能很好的解决关于颗粒的后处理问题了。一种是要想知道某个颗粒归属哪个网格就一直追踪这个颗粒的坐标，然后再确定网格；另一种是要知道某个网格内有哪几个粒子就像你说的遍历一遍网格确定粒子，会很慢。但可能还会遇到下面的问题，我在做CFD-DEM的时候在paraview中选取的粒子发现有两个编号，一个是：ID，一个是id，而且在多核计算的时候发现id会在跨越block时进行重新编号（后处理显示时），DPMFoam中应该不存在。

李东岳

@anubis 我当时也想过每个网格建立数组。在我这个算例里面（几千个网格几百个粒子）m个数组倒是可以接受。遍历粒子的时候，可以把网格数组的大小固定。但这一步骤，只是固定了数组大小。比如网格1的数组有3个，网格2的数组有5个，网格3的数组是2个。后续还要遍历一遍进行填充？

要是能遍历一次就好了

李东岳

@fubianhanshu 嗯，最慢的就是先来个遍历网格，然后每个网格里面遍历粒子这个循环套起来就很大了

anubis

@李东岳
遍历一次就可以了吧，我原来的想法是每次往数组最后添加元素（粒子的编号），不过刚搜了一下C++的数组大小是固定的（C++学的不太扎实），换成vector可能可以？

李东岳

@anubis vector只能加3个。如果数组固定的话，就需要一个比较大的数组，比如20 30个（针对我这个算例）。如果用不固定的数组，第一遍就需要确定数据元素数，然后第二遍添加元素。需要两次。如果一次能搞定就好了

gyzhangqm

粒子数量n较少，网格数量m较多时，可以对网格单元建立ADT树，然后确定每个粒子位于哪一个网格，复杂度nlog(m)

OItoCFD

@anubis 那就用位压缩，只需要建立一个m大小的一维数组。每个位置只存储一个该网格对应的粒子有无信息。比如一共十个粒子，网格23对应的状态就是0000100011，转换成十进制为35，只需要把35存进数组位置，因为一个整数35就代表了1，2，6三个位置处状态为1也就是这三个编号粒子存在。下次需要该网格23的粒子编号信息时候，取出35，再写个小程序转换进制判断一下，时间复杂度很低的，一般就是一点简单的位运算，就能知道换成二进制后第1，2，6个位置处是1.

李东岳

@gyzhangqm 我查了查ADT树这东西。叫二叉树？朋友圈里面有2个人说用二叉树。但我还不知道这玩意咋玩我研究研究

Jacobian

把粒子按x坐标排序，i是序号，形成xorder(i,1)存x坐标，xorder(i,2)存粒子编号；
同理，把粒子按y、z坐标排序，形成yorder(i,1..2)，zorder(i,1..2)
计算网格顶点的xmin，xmax，ymin，ymax，zmin，zmax

loop:循环网格{
      
       找到xmin，xmax在xorder(i,1)中的位置
       同理，处理ymin，ymax，zmin，zmax
       new 一个数组 possible(i)
     
       loop: 处理上面缩小的范围{
              把以上范围内的粒子编号xorder(i,2)，yorder(i,2),
              zorder(i,2)插入排 序到possible(i)
       }

       loop: 确定真正在网格中的粒子{
             possible(i)中有连续三个相同编号，判断是否真的在网格内,
             如果真在，就存起来
       }
}

至于每个网格里的粒子编号怎么存，搞一个数组叫data(i)一个接一个地存上面伪代码得到的编号，另一个数组叫position(j)存单元j在data(i)中的起始位置
不知道问题理解对不对，抛个砖试试看

马乔

用个List<dynamicList<label>> parcelsToCell, 然后遍历粒子，得到单元编号，然后将粒子塞到对应的dynamicList中。但是有空的。
或者用个std::Multimap<label>，cell编号就是key，value是粒子id。

同学博

       map<int,int> A;
       // particle/cell
      A.insert(pair<int,int>(1,1));
      for (map<int,int>::iterator iter = A.begin(); iter!=A.end(); ++iter)
      {
          if(iter->second==cellIDYouWant)
          {
              cout << iter->first << endl;
          }
      }

马乔

@同学博 map用迭代器访问效率并不高吧，而且最好不要直接访问value的吧，如果用value匹配，这样复杂度又上来了

夏雨天

对粒子并行，每个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是个好出发点

李东岳

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

网格遍历
粒子遍历
如果粒子恰好在网格内，就把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

同学博

@马乔