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

ir77

李老师及各位大佬好：.

我想要用自适应网格加密气液界面。

但不想按照常规AMR做法，仅加密0.001<α<0.999的网格，因为这样仅能加密气液界面的几层网格，我想多加密一些气液界面的网格。我考虑：能否定义一个标量场（例如该场名为regionforrefine），标记气液界面附近一定距离的网格，将这些网格的regionforrefine赋值为1，其他网格regionforrefine赋值为0，然后在dynamicMeshdict里指定加密0.5<regionforrefine<1.5的网格。我想这样就能多加密一些气液界面的网格（如图）

目前小弟在求解器的C文件里写了下面这一些代码，仅仅能够实现将气液界面单独标记出来，但如何将标记区域向内和向外扩张，我没有一点思路。

        if (alpha[i]>0&&alpha[i]<1)
           regionforrefine.field()[i]=1;
        else
           regionforrefine.field()[i]=0;

希望能够得到李老师和各位大佬的点拨

bestucan

把这四行代码写成个函数 forRefineFieldUpdate()

对你的代码搜 alpha，
碰见 alpha[i] = bala bala 这样的，就在后面放个 update 函数

就是：哪更新了 alpha的值，就在哪 update 一下。

或者：在 regionforrefine.field()[i]被调用之前更新

ir77

@bestucan bestucan大佬您好，我的那几行代码，仅能实现对气液界面单层网格的标记（如下图的红色网格，代表标记值1），但是我想要标记气液界面附近的多层网格：比如界面的内部5层外部5层的网格我都想把它标记下来。大佬您有相关建议吗

李东岳

@ir77 在标记气液界面附近一定距离的网格用于自适应网格加密中说：

比如界面的内部5层外部5层的网格我都想把它标记下来

即使在数学上，这个也很难定义。你怎么定义所谓的界面的内部5层外部5层的网格呢？

ir77

@李东岳
李老师好！
对于我的问题，我有两个思路，但限于水平太菜，两个思路都弄不出来：

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

有一个简单粗暴的方法。

原始单层的网格遍历两次：
第一次，给原始单层网格每个网格的左右+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楼的代码完成了代码修改，实现了对气液界面附近的更多网格的标记。
82bce252-6c29-43b8-aadf-80658d9cf8fd-QQ图片20220424200458.png
代码放在下面，供大家参考，我是初学者，代码写的不规范，大家见谅

      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]才可以正常使用。
大佬您能从经验上帮我判断一下可能是哪里出错了吗
27f57d6b-6b33-4fc0-b061-74860f37277a-K4_)8AHAPK1N9B4DZ}59YNM.png
$b8691bcd-39b0-46d7-8a88-46a41411da64-}93K@JY~XZT%[R0_]{$5V%6.png$
bac9c337-2063-4585-8d8d-fface8c9a54b-C%6S~[G}}R]Z274E7Y0K0`D.png

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的接口读取加密时间步间隔和加密范围，测试中发现程序只会在刚开始运行时读取这些常量。请问各位大佬，有没有什么办法能否让程序运行中也能读取这些常量呢？因为假设运行过程中我想改动这两个常量，目前只能把计算停下来改了后重新继续算，不能边算边改。