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的分享 -
@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是自定义的类,主要是函数的实现。
-
应该是精度的问题,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<>这个类的实现