标记气液界面附近一定距离的网格用于自适应网格加密

[bestucan]

有一个简单粗暴的方法。

原始单层的网格遍历两次：
第一次，给原始单层网格每个网格的左右+5的相邻网格都标记上。
第二次，给原始单层网格每个网格的上下+5的相邻网格都标记上。

为了第二次遍历时；还能找到“原始单层网格”，而不包括第一次遍历标记上的+5的网格；需要给“原始单层网格”再加个特殊标记。

这两次遍历会有重复标记的，不影响。

至于怎么找左右相邻的网格。。。应该能找到吧，这种那么规整的

[ir77]

bestucan大佬好，我能明白您的意思，脑子里也清楚整个流程，但是这一步*（给原始单层网格每个网格的左右+5的相邻网格都标记上）*我确实不知道该怎么用OpenFOAM语法表示出来

我能想到的只有我在第一楼最后提到的遍历网格，将0.001<α<0.999的网格赋值为1，但怎么将这个网格上下左右前后五个网格范围内的网格值也赋值为1，我不清楚该怎么弄

希望bestucan大佬能再花一点点时间帮我考虑一下这个问题

[bestucan]

@ir77 拿这个加工应该就行了

[ir77]

@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大佬费心了，谢谢大佬

[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，然后在第二步循环这个数组。

[ir77]

@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;

[ir77]

@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]才可以正常使用。
大佬您能从经验上帮我判断一下可能是哪里出错了吗

[bestucan]

可能你想要用的是 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 网页高亮变量名，看出现和变动的位置。

[ir77]

@bestucan 太谢谢大佬了！我研究研究

[ir77]

大家好，我来更新这个问题的相关进展。最近两天我查了些资料，基本上解决了这个问题，感谢大家的帮助，尤其是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的接口读取 加密时间步间隔 和 加密范围，测试中发现程序只会在刚开始运行时读取这些常量。请问各位大佬，有没有什么办法能否让程序运行中也能读取这些常量呢？因为假设运行过程中我想改动这两个常量，目前只能把计算停下来改了后重新继续算，不能边算边改。