粒子与网格归属问题

李东岳

@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()));
}

都是我云的

CFD中文网

粒子与网格归属问题