CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    OpenFOAM编程findCell的诡异问题,对同一坐标寻找cell结果不一致

    OpenFOAM
    5
    15
    1403
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • O
      OItoCFD 最后由 编辑

      各位老师同学好,最近写的一个代码的一个小部分用到了findcell去找点所在的cell,但是发现了一个奇怪的问题。对于两个点,我单独findcell出来结果是对的,但是在循环里面算出来这两个点后再findcell,结果不对。首先,代码我已经把无用的都删了后如下:

      #include "fvCFD.H"
      #include "meshSearch.H"
      
      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      
      int main(int argc, char *argv[])
      {
          argList::addBoolOption
          (
              "noOverwrite",
              "increment time to avoid overwriting initial fields"
          );
          #include "setRootCase.H"
          #include "createTime.H"
      
          // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      
          #include "createMesh.H"
      IOdictionary blockMeshDict
      (
          IOobject
          (
              "blockMeshDict",
              runTime.system(),
              mesh,
              IOobject::MUST_READ,
              IOobject::NO_WRITE
          )
      );  
      IOdictionary NetProperties
      (
          IOobject
          (
              "NetProperties",
              runTime.constant(),
              mesh,
              IOobject::MUST_READ,
              IOobject::NO_WRITE
          )
      );
      const dictionary& NetDict(NetProperties.subDict("deltaFunction"));
      const label  nx(blockMeshDict.lookup<scalar>("nx"));
      const label  ny(blockMeshDict.lookup<scalar>("ny"));
      const label  nz(blockMeshDict.lookup<scalar>("nz"));
      const scalar lx(blockMeshDict.lookup<scalar>("lx"));
      const scalar ly(blockMeshDict.lookup<scalar>("ly"));
      const scalar lz(blockMeshDict.lookup<scalar>("lz"));
      const scalar dx = lx/nx;
      const scalar dy = ly/ny;
      const scalar dz = lz/nz;
      const List<label> offSet_
      (
          NetDict.lookup("a")
      );
      
      meshSearch searchEngine_(mesh);
      
      vector tmpKnotPoint(1.5,0.72,1.359);
      vector point1(1.47727,0.694583,1.30181);
      vector point2(1.5,0.694583,1.30181);
      Pout << "point1 (1.47727 0.694583 1.30181) should be inside cell " << searchEngine_.findCell(point1) << endl;
      Pout << "point2(1.5 0.694583 1.30181) should be inside cell " << searchEngine_.findCell(point2) << endl;
      
      vector tmpOffSetPoint;
      label num = 0;
      forAll(offSet_, zi)
      {
          tmpOffSetPoint.z() = tmpKnotPoint.z() + offSet_[zi]*dz;
          forAll(offSet_, yi)
          {
              tmpOffSetPoint.y() = tmpKnotPoint.y() + offSet_[yi]*dy; 
              forAll(offSet_, xi)
              {
                  tmpOffSetPoint.x() = tmpKnotPoint.x() + offSet_[xi]*dx;
                  label cellId = searchEngine_.findCell(tmpOffSetPoint);
                  Pout << num << "  " << cellId << "  " << tmpOffSetPoint  << endl;
              }                    
          }
      }
      
          return 0;
      }
      

      输出如下:

      point1 (1.47727 0.694583 1.30181) should be inside cell 423586
      point2(1.5 0.694583 1.30181) should be inside cell 423587
      0  423587  (1.47727 0.694583 1.30181)
      0  423587  (1.5 0.694583 1.30181)
      0  423785  (1.47727 0.72 1.30181)
      0  423785  (1.5 0.72 1.30181)
      0  442595  (1.47727 0.694583 1.359)
      0  442595  (1.5 0.694583 1.359)
      0  442793  (1.47727 0.72 1.359)
      0  442793  (1.5 0.72 1.359)
      

      循环里输出的Cell 的Id对于这两个点(1.47727 0.694583 1.30181)和(1.5 0.694583 1.30181)都是423587,但是上面的输出直接输出这两个坐标findcell的结果,那就是一个423586,另一个是423587.

      难道是findcell用的方法不对,还是什么其他原因呢?请各位同学老师帮下忙。

      1 条回复 最后回复 回复 引用
      • O
        OItoCFD 最后由 编辑

        我把代码和测试算例精简后打包发在了网盘,如果有兴趣的老师和朋友可以帮忙看一下,谢谢!blockMesh后运行
        链接:https://pan.baidu.com/s/1HEy6w8MgbeXu9EsMcayung
        提取码:1111
        --来自百度网盘超级会员V4的分享

        1 条回复 最后回复 回复 引用
        • 李东岳
          李东岳 管理员 最后由 编辑

          代码上看起来是没问题。

          label cellId = searchEngine_.findCell(tmpOffSetPoint);

          这句话之前赋值为0然后再赋值呢?

          CFD高性能服务器 http://dyfluid.com/servers.html

          O 1 条回复 最后回复 回复 引用
          • 李东岳
            李东岳 管理员 最后由 编辑

            另外你没用过primitiveMesh的findCell()试试么,这种

            point a(1,2,3);
            mesh.findCell(a);
            

            我几年前记得在某个SCI里面看到过说meshSearch好像效率低,不过具体细节我忘了

            CFD高性能服务器 http://dyfluid.com/servers.html

            O 1 条回复 最后回复 回复 引用
            • O
              OItoCFD @李东岳 最后由 编辑

              @李东岳 在 OpenFOAM编程findCell的诡异问题,对同一坐标寻找cell结果不一致 中说:

              代码上看起来是没问题。

              label cellId = searchEngine_.findCell(tmpOffSetPoint);

              这句话之前赋值为0然后再赋值呢?

              老师,这样做也是一样的结果,还是重复的CellID

              1 条回复 最后回复 回复 引用
              • O
                OItoCFD @李东岳 最后由 编辑

                @李东岳 在 OpenFOAM编程findCell的诡异问题,对同一坐标寻找cell结果不一致 中说:

                另外你没用过primitiveMesh的findCell()试试么,这种

                point a(1,2,3);
                mesh.findCell(a);
                

                我几年前记得在某个SCI里面看到过说meshSearch好像效率低,不过具体细节我忘了

                老师,mesh search和mesh的findCell我都是用了的,一样的结果,都是循环外单独找那两个点是对的,循环内就不对,是重复的CellID

                马乔 1 条回复 最后回复 回复 引用
                • 马乔
                  马乔 副教授 @OItoCFD 最后由 编辑

                  @oitocfd 可能是精度问题?手动写入的有截断误差?刚好这个点位于face上呢?还有不管是primitiveMesh和meshSearch的findcell效率都挺低的,都要去遍历所有的cell,不知道写什么需求要这么写?

                  装逼没输过,吵架没赢过!

                  O 1 条回复 最后回复 回复 引用
                  • O
                    OItoCFD @马乔 最后由 编辑

                    @马乔 您好,这么讲的话,也许是可能是截断误差影响。然后我之所以用findCell,是因为我用到一个delta function,做一个点的物理量分布到周围网格范围内流体网格的distribution的操作,和点周围流体网格物理量插值到点的interpolation操作。就是grid to point and point to grid。需要找到点所在cell及其周围一圈两倍网格间距范围内的cell,不知道想实现这种插值分布函数,还有没有其他更简单的方法吗?请问您有什么看法吗?谢谢!

                    马乔 上级 2 条回复 最后回复 回复 引用
                    • 马乔
                      马乔 副教授 @OItoCFD 最后由 编辑

                      @oitocfd 我错了,meshSearch的查找效率很高,是primitiveMesh不行。你这是在写pic吗?

                      装逼没输过,吵架没赢过!

                      O 1 条回复 最后回复 回复 引用
                      • 上级
                        上级 @OItoCFD 最后由 编辑

                        @oitocfd 之前写过一个小函数,具体是在dpmfoam里用到的,是寻找一个颗粒所在的cell label以及这个cell周围的一圈cell label,不知道能不能帮到你。

                        template<class Type>
                        const Foam::labelList& Foam::interpolationCellAround<Type>::cellpointCells
                        (
                            const label celli,
                            DynamicList<label>& storage
                        ) const
                        {
                                const labelList& cPoints = this->psi_.mesh().cellPoints()[celli];
                        
                                storage.clear();
                        
                                forAll(cPoints, i)
                                {
                                    const label pointi = cPoints[i];
                        
                                    const labelList& pointcellsi = this->psi_.mesh().pointCells()[pointi];
                        
                                    forAll(pointcellsi, j)
                                    {
                                        storage.append(pointcellsi[j]);
                                    }
                                }
                        
                                // Filter duplicates
                                if (storage.size() > 1)
                                {
                                    sort(storage);
                        
                                    label n = 1;
                                    for (label i = 1; i < storage.size(); i++)
                                    {
                                        if (storage[i-1] != storage[i])
                                        {
                                            storage[n++] = storage[i];
                                        }
                                    }
                        
                                    // truncate addressed list
                                    storage.setSize(n);
                                }
                        
                                return storage;
                        }
                        
                        template<class Type>
                        const Foam::labelList& Foam::interpolationCellAround<Type>::cellpointCells(const label celli) const
                        {
                            return cellpointCells(celli, labels_);
                        }
                        
                        

                        interpolationCellAround是自定义的类,主要是函数的实现。

                        O 1 条回复 最后回复 回复 引用
                        • O
                          OItoCFD @马乔 最后由 编辑

                          @马乔 您好,算是吧 就是每个时间步 我要知道我每个拉格朗日点在哪个cell,然后把力按照cell距离的远近按照delta函数分布到周围所有cell里 我的固体是拉格朗日点代表的 固体可变形 所以每个时间步都要所有点来一次findCell

                          1 条回复 最后回复 回复 引用
                          • O
                            OItoCFD @上级 最后由 编辑

                            @上级 在 OpenFOAM编程findCell的诡异问题,对同一坐标寻找cell结果不一致 中说:

                            mesh().pointCells()

                            好的 我学习一下 谢谢

                            马乔 1 条回复 最后回复 回复 引用
                            • 马乔
                              马乔 副教授 @OItoCFD 最后由 编辑

                              @oitocfd 如果粒子追踪的话可以尝试下lagrangian库,当前网格周围一圈网格可以直接用primitiveMesh.cellCells()

                              装逼没输过,吵架没赢过!

                              1 条回复 最后回复 回复 引用
                              • Tong
                                Tong 最后由 编辑

                                应该是精度的问题,polyMesh 的findCell 的原型为:findCell(const point & p,const cellDecomposition decompMode = CELL_TETS) 默认是使用polyMesh里的八叉树进行空间搜索的,vector每个点是一个双精度值,有效数字大约有11位置,而Info只输出了部分有效数字,感觉你直接这样写出来,就会损失精度

                                vector tmpKnotPoint(1.5,0.72,1.359);
                                vector point1(1.47727,0.694583,1.30181);
                                vector point2(1.5,0.694583,1.30181);
                                

                                如果网格尺寸很小的话可能造成差异,
                                还有就是,我的记得polyMesh 网格索引树默认的精度是1e-3 如果你的网格尺寸小的话你可以调整一下,具体的你可以去看indexedOctree<>这个类的实现

                                O 1 条回复 最后回复 回复 引用
                                • O
                                  OItoCFD @Tong 最后由 编辑

                                  @tong 好的好的 感谢感谢

                                  1 条回复 最后回复 回复 引用
                                  • First post
                                    Last post