Skip to content
  • 最新
  • 版块
  • 东岳流体
  • 随机看[请狂点我]
皮肤
  • Light
  • Cerulean
  • Cosmo
  • Flatly
  • Journal
  • Litera
  • Lumen
  • Lux
  • Materia
  • Minty
  • Morph
  • Pulse
  • Sandstone
  • Simplex
  • Sketchy
  • Spacelab
  • United
  • Yeti
  • Zephyr
  • Dark
  • Cyborg
  • Darkly
  • Quartz
  • Slate
  • Solar
  • Superhero
  • Vapor

  • 默认(不使用皮肤)
  • 不使用皮肤
折叠
CFD中文网

CFD中文网

  1. CFD中文网
  2. OpenFOAM
  3. OpenFOAM编程findCell的诡异问题,对同一坐标寻找cell结果不一致

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

已定时 已固定 已锁定 已移动 OpenFOAM
15 帖子 5 发布者 9.7k 浏览
  • 从旧到新
  • 从新到旧
  • 最多赞同
回复
  • 在新帖中回复
登录后回复
此主题已被删除。只有拥有主题管理权限的用户可以查看。
  • O 离线
    O 离线
    OItoCFD
    写于2021年11月22日 21:51 最后由 编辑
    #1

    各位老师同学好,最近写的一个代码的一个小部分用到了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 离线
    O 离线
    OItoCFD
    写于2021年11月22日 21:57 最后由 编辑
    #2

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

    1 条回复 最后回复
  • 李 在线
    李 在线
    李东岳 管理员
    写于2021年11月22日 23:49 最后由 编辑
    #3

    代码上看起来是没问题。

    label cellId = searchEngine_.findCell(tmpOffSetPoint);

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

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    O 1 条回复 最后回复 2021年11月23日 10:04
  • 李 在线
    李 在线
    李东岳 管理员
    写于2021年11月23日 01:19 最后由 编辑
    #4

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

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

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

    http://dyfluid.com/index.html
    需要帮助debug算例的看这个 https://cfd-china.com/topic/8018

    O 1 条回复 最后回复 2021年11月23日 10:09
  • O 离线
    O 离线
    OItoCFD
    在 2021年11月23日 10:04 中回复了 李东岳 最后由 编辑
    #5

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

    代码上看起来是没问题。

    label cellId = searchEngine_.findCell(tmpOffSetPoint);

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

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

    1 条回复 最后回复
  • O 离线
    O 离线
    OItoCFD
    在 2021年11月23日 10:09 中回复了 李东岳 最后由 编辑
    #6

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

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

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

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

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

    马乔马 1 条回复 最后回复 2021年11月24日 04:00
  • 马乔马 离线
    马乔马 离线
    马乔 大神
    在 2021年11月24日 04:00 中回复了 OItoCFD 最后由 编辑
    #7

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

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

    O 1 条回复 最后回复 2021年11月24日 09:24
  • O 离线
    O 离线
    OItoCFD
    在 2021年11月24日 09:24 中回复了 马乔 最后由 编辑
    #8

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

    马乔马 上级上 2 条回复 最后回复 2021年12月14日 02:10
  • 马乔马 离线
    马乔马 离线
    马乔 大神
    在 2021年12月14日 02:10 中回复了 OItoCFD 最后由 编辑
    #9

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

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

    O 1 条回复 最后回复 2021年12月14日 09:45
  • 上级上 离线
    上级上 离线
    上级
    在 2021年12月14日 05:58 中回复了 OItoCFD 最后由 编辑
    #10

    @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 条回复 最后回复 2021年12月14日 09:56
  • O 离线
    O 离线
    OItoCFD
    在 2021年12月14日 09:45 中回复了 马乔 最后由 编辑
    #11

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

    1 条回复 最后回复
  • O 离线
    O 离线
    OItoCFD
    在 2021年12月14日 09:56 中回复了 上级 最后由 编辑
    #12

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

    mesh().pointCells()

    好的 我学习一下 谢谢

    马乔马 1 条回复 最后回复 2021年12月24日 01:51
  • 马乔马 离线
    马乔马 离线
    马乔 大神
    在 2021年12月24日 01:51 中回复了 OItoCFD 最后由 编辑
    #13

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

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

    1 条回复 最后回复
  • TongT 离线
    TongT 离线
    Tong
    写于2021年12月30日 11:17 最后由 编辑
    #14

    应该是精度的问题,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 条回复 最后回复 2022年1月1日 12:46
  • O 离线
    O 离线
    OItoCFD
    在 2022年1月1日 12:46 中回复了 Tong 最后由 编辑
    #15

    @tong 好的好的 感谢感谢

    1 条回复 最后回复
2021年11月22日 21:51

1/15

2021年11月22日 21:51

未读 14
2022年1月1日 12:46
  • 登录

  • 登录或注册以进行搜索。
1 / 15
  • 第一个帖子
    1/15
    最后一个帖子
0
  • 最新
  • 版块
  • 东岳流体
  • 随机看[请狂点我]