OpenFOAM编程findCell的诡异问题，对同一坐标寻找cell结果不一致
-
各位老师同学好，最近写的一个代码的一个小部分用到了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用的方法不对，还是什么其他原因呢？请各位同学老师帮下忙。
-
我把代码和测试算例精简后打包发在了网盘，如果有兴趣的老师和朋友可以帮忙看一下，谢谢！blockMesh后运行
链接：https://pan.baidu.com/s/1HEy6w8MgbeXu9EsMcayung
提取码：1111
--来自百度网盘超级会员V4的分享 -
代码上看起来是没问题。
label cellId = searchEngine_.findCell(tmpOffSetPoint);
这句话之前赋值为0然后再赋值呢？
-
另外你没用过primitiveMesh的findCell()试试么，这种
point a(1,2,3); mesh.findCell(a);我几年前记得在某个SCI里面看到过说meshSearch好像效率低，不过具体细节我忘了
-
@李东岳 在 OpenFOAM编程findCell的诡异问题，对同一坐标寻找cell结果不一致 中说：
代码上看起来是没问题。
label cellId = searchEngine_.findCell(tmpOffSetPoint);
这句话之前赋值为0然后再赋值呢？
老师，这样做也是一样的结果，还是重复的CellID
-
@李东岳 在 OpenFOAM编程findCell的诡异问题，对同一坐标寻找cell结果不一致 中说：
另外你没用过primitiveMesh的findCell()试试么，这种
point a(1,2,3); mesh.findCell(a);我几年前记得在某个SCI里面看到过说meshSearch好像效率低，不过具体细节我忘了
老师，mesh search和mesh的findCell我都是用了的，一样的结果，都是循环外单独找那两个点是对的，循环内就不对，是重复的CellID
-
@oitocfd 可能是精度问题？手动写入的有截断误差？刚好这个点位于face上呢？还有不管是primitiveMesh和meshSearch的findcell效率都挺低的，都要去遍历所有的cell，不知道写什么需求要这么写？
-
@马乔 您好，这么讲的话，也许是可能是截断误差影响。然后我之所以用findCell，是因为我用到一个delta function，做一个点的物理量分布到周围网格范围内流体网格的distribution的操作，和点周围流体网格物理量插值到点的interpolation操作。就是grid to point and point to grid。需要找到点所在cell及其周围一圈两倍网格间距范围内的cell，不知道想实现这种插值分布函数，还有没有其他更简单的方法吗？请问您有什么看法吗？谢谢！
-
@oitocfd 我错了，meshSearch的查找效率很高，是primitiveMesh不行。你这是在写pic吗？
-
@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是自定义的类，主要是函数的实现。
-
@马乔 您好，算是吧 就是每个时间步 我要知道我每个拉格朗日点在哪个cell，然后把力按照cell距离的远近按照delta函数分布到周围所有cell里 我的固体是拉格朗日点代表的 固体可变形 所以每个时间步都要所有点来一次findCell
-
-
@oitocfd 如果粒子追踪的话可以尝试下lagrangian库，当前网格周围一圈网格可以直接用primitiveMesh.cellCells()
-
应该是精度的问题，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<>这个类的实现 -
@tong 好的好的 感谢感谢
