标记气液界面附近一定距离的网格用于自适应网格加密
-
-
@李东岳
李老师好！
对于我的问题，我有两个思路，但限于水平太菜，两个思路都弄不出来：1.标记气液界面附近的一定层数的网格，也就是前面说的“内部5层外部5层的网格”。我找到的比较相关的代码位于AMR源文件的dynamicRefineFvMesh.C里：
// Extend with a buffer layer to prevent neighbouring points // being unrefined. for (label i = 0; i < nBufferLayers; ++i) { extendMarkedCells(refineCell); }这里的nBufferLayers是执行AMR时在dynamicMeshDict里填写的数字，可以通过这段代码控制缓冲层的层数。我想的是将气液界面标记成一个标量场（存在气液界面的网格赋值为1，其他网格赋值为0），并借鉴这段代码将标记区向内和向外“延拓”几个网格。但水平太菜，尝试了多次未能成功。
2.标记气液界面附近的的一定距离的网格。我想的是定义一个标量场，标量场的初始值都为0，然后对网格做遍历，只要网格到气液界面（0.001<α<0.999）的距离小于一定距离则该网格标量赋值为1。但是这种遍历很复杂我不知道该怎么写。
我的水平很菜，如果理解有偏差，表述有不清楚的地方望李老师指出
谢谢李老师。 -
-
@bestucan 多谢大佬费心，我昨天也看到了这个帖子，昨天我参照这个帖子做出了一个初级版本，能够实现对气液界面附近三层网格的标记。
我采用的代码是这样的：const labelList& neighbour = mesh.cellCells()[i]; scalar sum(0); forAll(neighbour, cellI) { sum += alpha1[neighbour[cellI]]; } sum /= neighbour.size(); if (sum>0.0001&&sum<0.9999) regionforinterface.field()[i]=1; else regionforinterface.field()[i]=0;后来，我想给它添加一个嵌套进去，那样就可以加密气液界面附近五层网格了，采用如下代码：
const labelList& neighbour = mesh.cellCells()[i]; for(int j=0; j<neighbour.size(); j++) { const labelList& neighbour2= mesh.cellCells()[neighbour[j]]; } scalar sum(0); forAll(neighbour2, cellI) { sum += alpha1[neighbour2[cellI]]; } sum /= neighbour2.size(); if (sum>0.0001&&sum<0.9999) regionforinterface.field()[i]=1; else regionforinterface.field()[i]=0;但是编译的时候提示neighbour2未申明：
error: ‘neighbour2’ was not declared in this scope forAll(neighbour2, cellI)
bestucan大佬，能麻烦您帮我看下第二段代码吗，我不知道该怎么改
另外，bestucan大佬，我的第一段代码是参照东岳老师之前的一个帖子改的，是通过对目标网格所有相邻网格的alpha场做平均，如果平均值位于标记范围内则该目标网格标记为1。
我想，能否改一下代码，使得只要目标网格的任意一个相邻网格alpha位于标记范围内则该目标网格标记为1。感觉应该用if判断语句，但这里的判断对象又是多个网格中的alpha值，我不知道该怎么写这段代码，麻烦大佬指点一下迷津。
bestucan大佬费心了，谢谢大佬 -
@ir77 因为你的
const labelList& neighbour2声明在 for 循环里，所以出了 for 循环无法使用挪外面，和 neighbour 放一起声明，然后在 for 里赋值
const labelList & neighbour = mesh.cellCells()[i]; const labelList & neighbour2; for (int j = 0; j < neighbour.size(); j++) { neighbour2 = mesh.cellCells()[neighbour[j]]; }得两步：
第一步，找出目标网格。
第二部，去判断目标网格并赋值找出目标网格这一步，新建一个数组，这个数组里放目标网格的网格ID，然后在第二步循环这个数组。
-
@bestucan bestucan 大佬您好，我提到的
只要目标网格的任意一个相邻网格alpha位于标记范围内则该目标网格标记为1
我通过下面这段代码解决了
const labelList& neighbour = mesh.cellCells()[i]; scalar sumforinterface(0); forAll(neighbour, cellI) { if (alpha1[neighbour[cellI]]>0.0001&&alpha1[neighbour[cellI]]<0.9999) sumforinterface += 1; else sumforinterface += 0; } if (sumforinterface>0) regionforinterface.field()[i]=1; else regionforinterface.field()[i]=0;然后我按照您建议的
因为你的 const labelList& neighbour2声明在 for 循环里，所以出了 for 循环无法使用
挪外面，和 neighbour 放一起声明，然后在 for 里赋值将const labelList & neighbour2;放到了外边，但是报错error: ‘neighbour2’ declared as reference but not initialized
好像是说需要初始化，我不知道该怎么弄
后来我在cfdonline上找到了一条和我问题类似的帖子。参照18楼的代码完成了代码修改，实现了对气液界面附近的更多网格的标记。
代码放在下面，供大家参考，我是初学者，代码写的不规范，大家见谅
const labelList& neighbour = mesh.cellCells()[i]; labelHashSet setNBCells(1); for(int j=0; j<neighbour.size(); j++) { const labelList & neighbour2= mesh.cellCells()[neighbour[j]]; setNBCells.insert(neighbour2); } labelList neighbour2=setNBCells.toc(); scalar sumforinterface(0); forAll(neighbour2, cellI) { if (alpha1[neighbour2[cellI]]>0.0001&&alpha1[neighbour2[cellI]]<0.9999) sumforinterface += 1; else sumforinterface += 0; } if (sumforinterface>0) regionforinterface.field()[i]=1; else regionforinterface.field()[i]=0; -
@bestucan 对了大佬，我还碰到了一个小小的问题：
@ir77 在 标记气液界面附近一定距离的网格用于自适应网格加密 中说：
const labelList& neighbour = mesh.cellCells()[i];
scalar sumforinterface(0);
forAll(neighbour, cellI)
{
if (alpha1[neighbour[cellI]]>0.0001&&alpha1[neighbour[cellI]]<0.9999)
sumforinterface += 1;
else
sumforinterface += 0;
}
if (sumforinterface>0)
regionforinterface.field()[i]=1;
else
regionforinterface.field()[i]=0;我把这段代码第一句中的mesh.cellCells()[i]分别改成了mesh.pointCells()[i]和mesh.edgeCells()[i]，运行结果却非常糟糕（如下图），似乎只有mesh.cellCells()[i]才可以正常使用。
大佬您能从经验上帮我判断一下可能是哪里出错了吗
![b8691bcd-39b0-46d7-8a88-46a41411da64-}93K@JY~XZT%[R0_]{$5V%6.png](/assets/uploads/files/1650803378336-b8691bcd-39b0-46d7-8a88-46a41411da64-93k-jy-xzt-r0_-5v-6.png)
![bac9c337-2063-4585-8d8d-fface8c9a54b-C%6S~[G}}R]Z274E7Y0K0`D.png](/assets/uploads/files/1650803395252-bac9c337-2063-4585-8d8d-fface8c9a54b-c-6s-g-r-z274e7y0k0-96-d.png)
-
可能你想要用的是 mesh.cellPoints() 和 mesh.cellEdges()？
用一圈网格得到点，再用点得到相邻网格？
用一圈网格得到边，再用边得到相邻网格？参考：https://jibranhaider.com/blog/mesh-information-in-openfoam/
如果要从代码里看：：
cellCells这个函数在这里：
[1] https://cpp.openfoam.org/v9/primitiveMeshCellCells_8C_source.html
第100行，它返回的指针 ccptr 是经过 calCellCells() 新建的。calCellCells()这个函数就在 [1] 的31行。73行新建，78行到95行寻找、赋值。edgeCells()这个函数在这里：
[2] https://cpp.openfoam.org/v9/primitiveMeshEdgeCells_8C_source.html
第32行，它返回的指针 ecptr 是经过 invertManyToMany() 更新的。而 invertManyToMany() 这个函数在
[3] https://cpp.openfoam.org/v9/ListOpsTemplates_8C_source.html
第 433 行。invertManyToMany() 中有两个 forAll 函数，函数中的 pointEdges 是函数的传入参数，也就是 [2] 中50行调用该函数时的第二个参数 cellEdges()。这个 cellEdges() 在
[4] https://cpp.openfoam.org/v9/primitiveMeshCellEdges_8C_source.html
中第118行。118行调用的函数 calcCellEdges() 就在 [4] 中32行。calcCellEdges() 中对 cellEdges() 返回指针 ceptr 的操作在104 和105行（104行用63行新建的ce初始化ceptr，105行用别名 cellEdgeAddr 对 ceptr 操作）。这个操作过程和[1]中73，74行类似。对比[4]中67，68和[1]中63，64。都有 own 和nei。不同的是[4]中对 own 和 nei 的两个循环，把结果都添加到 curCellEdges 了（83和99行），而 curCellEdges 是 ce 的引用，ce又在104行用于初始化ceptr。但是，它储存的方式和 [1] 中不一样。
有点晕了
大概就是这么看。用网页右上方搜索框搜函数，用 ctrl + F 网页高亮变量名，看出现和变动的位置。
-
大家好，我来更新这个问题的相关进展。最近两天我查了些资料，基本上解决了这个问题，感谢大家的帮助，尤其是bestucan大佬。
参考帖子关于相邻cell我写了个遍历进行暴力求解。
（1）代码思路：对所有网格做遍历，找到气液界面网格，针对每一个气液界面网格都对所有网格再做一次遍历，计算两网格之间的距离，对满足距离要求的网格在markRegionNearInterface这个场上+1。最后，可以利用markRegionNearInterface>=1来筛选出所有满足条件的网格。
（2）这段代码计算起来实在是费时间，因为我的需求(见贴名)不需要每一个时间步都运行这段代码，所以我在代码最外层套了一个时间步判断，可以实现多个时间步运行一次该代码。
（3）我从其他地方抄了点代码写了两个接口，分别从transportProperties读取 加密时间步间隔 和 加密范围，其中加密范围是依据气泡半径写的无量纲范围，所以需要配合前面代码(未在本段代码中给出)求出的气泡半径radius.value()定义标记区域。代码如下
//read timestep interval required for the following code, as the following code is very time consuming so I needn't it run for every timestep. const label calculateTimestepInterval = transportProperties.get<label>("calculateTimestepInterval"); //read information from transportProperties so that we can modify the size of marked region easily. scalar markRegionNearInterfaceSize = transportProperties.get<scalar>("markRegionNearInterfaceSize"); if (mesh.time().timeIndex() % calculateTimestepInterval == 0)//judge whether the timestep satisfy the timestep interval criterion. { markRegionNearInterface.field()=0;//reset the field to be zero first. Of course such field is claimed in createFields.H first. forAll(mesh.C(), i)//Cycle all meshes { if (alpha1[i]>0.0001&&alpha1[i]<0.9999)//judge whether the cell contains gas-liquid interface { vector centerOfSurfacePoint = mesh.C()[i];//store center information of such a interface cell forAll(mesh.C(), cellI)//Cycle all meshes again { scalar offset = mag(centerOfSurfacePoint - mesh.C()[cellI]);//distance between any cells and the interface cell if (offset <= markRegionNearInterfaceSize*radius.value())//judge the criterion, note the radius is calculated before with code which is not shown here { markRegionNearInterface.field()[cellI] += 1;//any cells which meet the criterion will be marked with value>0 } } } } }运行结果如下
BTW
（1） 因为我是个beginner，代码是东抄抄西抄抄弄出来的，上述代码可能有些冗余，如果大佬们有时间，想麻烦大佬们帮我看看我这些代码能从哪些角度优化优化。代码中的变量类型scalar vector label之类都是我一次次试错试出来的，也不知道这些变量类型用的合理不合理。
（2）我在代码中写了两个链接到transportProperties的接口读取 加密时间步间隔 和 加密范围，测试中发现程序只会在刚开始运行时读取这些常量。请问各位大佬，有没有什么办法能否让程序运行中也能读取这些常量呢？因为假设运行过程中我想改动这两个常量，目前只能把计算停下来改了后重新继续算，不能边算边改。
