OpenFOAM-v2006重叠网格挖洞问题研究以及overset代码解析（八）walkFront函数-卡核

在我之前的博客里最后提到一个解决方案是对walkFront函数进行改动的。一个叫louisgag的前辈在github上贴了他的代码。这次想要学习walkFront函数以及前辈的nPushWalkFront函数，为了弄明白这个函数在整一个重叠网格中担当的角色。这一函数比较重要，会尽量讲解每一行代码的意思（我理解的，如果有误欢迎指正！）

看到walkFront的参数walkFront(layerRelax, allStencil, allCellTypes, allWeight);其函数定义时参数：

其中，layerRelax是用户定义的标量，起作用的地方在后面的seedCell函数里。allStencil是上一篇博文提到的markDonors输出的下标列表类型变量。allCellType是在markPatchesAsHoles（在我前面五篇博文里有介绍）这个函数里输出的，在这个函数中还会进行进一步的改动。最后是allWeight，是一个标量场，其在调用walkFront前被定义：scalarField allWeight(mesh_.nCells(), Zero);此标量场是储存每一个网格的权重信息的，至于这个权重是用来干什么的，接着看下去。

    // Current front
    bitSet isFront(mesh_.nFaces());

    const fvBoundaryMesh& fvm = mesh_.boundary();

fvm调用了网格的边界mesh.boundary()；isFront则是用于标记挖洞之后的HOLE外面一圈的面。bitSet相当于一个只含有真假两种状态的列表。这里IsFront初始化时容量设置成所有面的数量，并初始化为0，对于需要的面会在后面标记为1。Front是一层需要被标记的插值层（buffer interpolation layer(s) around holes）。需要标记的区域有两类，一是OVERSET边界上的面，二是背景网格挖洞之后HOLE外面的区域。

1. OVERSET PATCH：

    // 'overset' patches

    forAll(fvm, patchI)
    {
        if (isA<oversetFvPatch>(fvm[patchI]))
        {
            const labelList& fc = fvm[patchI].faceCells();
            forAll(fc, i)
            {
                label cellI = fc[i];
                if (allCellTypes[cellI] == INTERPOLATED)
                {
                    // Note that acceptors might have been marked hole if
                    // there are no donors in which case we do not want to
                    // walk this out. This is an extreme situation.
                    isFront.set(fvm[patchI].start()+i);
                }
            }
        }
    }

如果边界上的某个patch为overset类型，那么这个patch上的面会被标记成Front。这个区域是两层网格的交界处，所以即使没有挖洞也必须标记为Front。这里还另外判断面上对应的cell是不是别标记为INTERPOLATED，是为了防止出现极端的情况：如果某一单元没有别哪个单元贡献的话，那么这个单元可能会被标记为HOLE，这种情况可能出现在重叠网格的外边界陷入到wall边界里面时。

2. outside of HOLE region：

    // Outside of 'hole' region
    {
        const labelList& own = mesh_.faceOwner();
        const labelList& nei = mesh_.faceNeighbour();

        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
        {
            label ownType = allCellTypes[own[faceI]];
            label neiType = allCellTypes[nei[faceI]];
            if
            (
                 (ownType == HOLE && neiType != HOLE)
              || (ownType != HOLE && neiType == HOLE)
            )
            {
                //Pout<< "Front at face:" << faceI
                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
                isFront.set(faceI);
            }
        }

        labelList nbrCellTypes;
        syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);

        for
        (
            label faceI = mesh_.nInternalFaces();
            faceI < mesh_.nFaces();
            faceI++
        )
        {
            label ownType = allCellTypes[own[faceI]];
            label neiType = nbrCellTypes[faceI-mesh_.nInternalFaces()];

            if
            (
                 (ownType == HOLE && neiType != HOLE)
              || (ownType != HOLE && neiType == HOLE)
            )
            {
                //Pout<< "Front at coupled face:" << faceI
                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
                isFront.set(faceI);
            }
        }
    }

这一块是标记背景网格里面HOLE外面一圈的Front的，其判断主要用：

if((ownType == HOLE && neiType != HOLE)
        || (ownType != HOLE && neiType == HOLE))

来实现。上部分是对背景网格内部面进行遍历的，而下部分则是对背景网格边界面进行遍历，总之每一次都要遍历背景网格的所有的面。这样，isFront就标记完毕，接下来是处理权重。

    // Current interpolation fraction
    scalarField fraction(mesh_.nFaces(), Zero);

    forAll(isFront, faceI)
    {
        if (isFront.test(faceI))
        {
            fraction[faceI] = 1.0;
        }
    }

这一段定义了一个标量场fraction，其大小和网格面元数量相等，且初始化为0，并根据isFront来设置某一个face上的值为1。这里用到的test是bool返回类型函数，返回faceI位置时的位值。在前面的代码中特定的面上isFront是1，所以这里这些面对应的fraction也会设置为1。

while (returnReduce(isFront.any(), orOp<bool>()))
    {
        // Interpolate cells on front
        bitSet newIsFront(mesh_.nFaces());
        scalarField newFraction(fraction);
        forAll(isFront, faceI)
        {
            if (isFront.test(faceI))
            {
                label own = mesh_.faceOwner()[faceI];
                if (allCellTypes[own] != HOLE)
                {
                    if (allWeight[own] < fraction[faceI])
                    {
                        // Cell wants to become interpolated (if sufficient
                        // stencil, otherwise becomes hole)
                        if (allStencil[own].size())
                        {
                            allWeight[own] = fraction[faceI];
                            allCellTypes[own] = INTERPOLATED;
                            // Add faces of cell (with lower weight) as new
                            // front
                            seedCell
                            (
                                own,
                                fraction[faceI]-layerRelax,
                                newIsFront,
                                newFraction
                            );
                        }
                        else
                        {
                            allWeight[own] = 0.0;
                            allCellTypes[own] = HOLE;
                            // Add faces of cell as new front
                            seedCell
                            (
                                own,
                                1.0,
                                newIsFront,
                                newFraction
                            );
                        }
                    }
                }
                if (mesh_.isInternalFace(faceI))
                {
                    label nei = mesh_.faceNeighbour()[faceI];
                    if (allCellTypes[nei] != HOLE)
                    {
                        if (allWeight[nei] < fraction[faceI])
                        {
                            if (allStencil[nei].size())
                            {
                                allWeight[nei] = fraction[faceI];
                                allCellTypes[nei] = INTERPOLATED;
                                seedCell
                                (
                                    nei,
                                    fraction[faceI]-layerRelax,
                                    newIsFront,
                                    newFraction
                                );
                            }
                            else
                            {
                                allWeight[nei] = 0.0;
                                allCellTypes[nei] = HOLE;
                                seedCell
                                (
                                    nei,
                                    1.0,
                                    newIsFront,
                                    newFraction
                                );
                            }
                        }
                    }
                }
            }
        }

        syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
        syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());

        isFront.transfer(newIsFront);
        fraction.transfer(newFraction);
    }

这是一个while的循环块：while (returnReduce(isFront.any(), orOp<bool>()))，isFront.any()是bool返回函数，当isFront里面只要有一个值是1时就会返回1；orOp<bool>()或运算，而returnReduce则是// Reduce using either linear or tree communication schedule ，实际上还是一个reduce函数。因为本人代码基础不佳，这里不能完全解读开来。但毕竟是一个while循环，我想当isFront所有的位都置0的时候，这一循环就终止了。为什么要有这个循环呢？我们先看里面的内容。它先是定义了两个新的量newIsFront和newFraction。可以推断出这个是用来替代原来的变量的。forAll(isFront, faceI)，对所有的面进行遍历，然后进行四层if判断，分别是：

1. 该面上对应的isFront的值是否为1？

2. 面的owner对应的cell的类型是不是非HOLE？或者该面是不是内部面？

3. allWeight是否小于fraction？

4. allStencil是否有大小？

newIsFront和newFraction会在if内的seedCell上进行改动。seedCell代码如下：

void Foam::cellCellStencils::inverseDistance::seedCell
(
    const label cellI,
    const scalar wantedFraction,
    bitSet& isFront,
    scalarField& fraction
) const
{
    const cell& cFaces = mesh_.cells()[cellI];
    forAll(cFaces, i)
    {
        label nbrFacei = cFaces[i];
        if (fraction[nbrFacei] < wantedFraction)
        {
            fraction[nbrFacei] = wantedFraction;
            isFront.set(nbrFacei);
        }
    }
}

结合函数调用时的参数，可以看到，最后对应面的newIsFront会被set，newFraction也会设为wantedFraction，这个值在调用时的第二个参数里。这个值被设置成fraction[faceI]-laterRelax或1。事实上，laterRelax一般默认为一个极小的值。在seedCell函数里，调用了这个体元的所有面，如果其中一个面的fraction低于wantedFraction，那么这个面的fraction就会等于wantedFraction，也就是说取二者较大的值，且此时newIsFront也会被set。

总的来看，这个walkFront函数的主角是面元face，而我们要进行标记的是cell。每一个face都会对应两个cell分别是own和nei，判断其是否要进行插值运算，则是根据allStencil是否有大小。而allStencil是在markDonors进行定义的。在markDonors里，如果两层网格里的两个单元相互有映射关系且其中任何一个都不是HOLE类型，那么就会有一个donor，也就是说allStencil就会有大小，在这里加上isFront的条件，使得这一个单元会被标记为INTERPOLATED类型，作为HOLE外层或者overset类型边界的插值单元，通过各种插值格式进行不同网格层之间的信息传输。

在循环体结尾还有四行代码：

syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());

isFront.transfer(newIsFront);
fraction.transfer(newFraction);

syncFaceList是 //- Synchronize values on all mesh faces.，然后transfer是把newIsFront替换掉isFront。至于为什么要这样做，以及为什么要进行循环，我想可能是和程序的优化有关系，这就不深究了。

挖洞的问题在各论坛里反响已久。其中louisgag提出了自己的解决方案：通过修改walkFront函数来使得挖洞可以往外继续进行。本篇博文介绍他的做法。

放上相关的代码：

int iPushFront = nPushFront;

    while (iPushFront > 0 && returnReduce(isFront.any(), orOp<bool>()))
    {
        // Interpolate cells on front
        bitSet newIsFront(mesh_.nFaces());
        //~ scalarField newFraction(fraction);
        forAll(isFront, faceI)
        {
            if (isFront.test(faceI))
            {
                label own = mesh_.faceOwner()[faceI];
                if (allCellTypes[own] != HOLE)
                {
                    //~ if (allWeight[own] < fraction[faceI])
                    //~ {
                        // Cell wants to become interpolated (if sufficient
                        // stencil, otherwise becomes hole)
                        //~ if (allStencil[own].size())
                        //~ {
                            //~ allWeight[own] = fraction[faceI];
                            allCellTypes[own] = HOLE;
                            newIsFront.set(mesh_.cells()[own]);
                            // Add faces of cell (with lower weight) as new
                            // front
                            //~ seedCell
                            //~ (
                                //~ own,
                                //~ fraction[faceI]-layerRelax,
                                //~ newIsFront,
                                //~ newFraction
                            //~ );
                        //~ }
                        //~ else
                        //~ {
                            //~ allWeight[own] = 0.0;
                            //~ allCellTypes[own] = HOLE;
                            //~ // Add faces of cell as new front
                            //~ seedCell
                            //~ (
                                //~ own,
                                //~ 1.0,
                                //~ newIsFront,
                                //~ newFraction
                            //~ );
                        //~ }
                    //~ }
                }
                if (mesh_.isInternalFace(faceI))
                {
                    label nei = mesh_.faceNeighbour()[faceI];
                    if (allCellTypes[nei] != HOLE)
                    {
                        //~ if (allWeight[nei] < fraction[faceI])
                        //~ {
                            //~ if (allStencil[nei].size())
                            //~ {
                                //~ allWeight[nei] = fraction[faceI];
                                allCellTypes[nei] = HOLE;
                                newIsFront.set(mesh_.cells()[nei]);
                                //~ seedCell
                                //~ (
                                    //~ nei,
                                    //~ fraction[faceI]-layerRelax,
                                    //~ newIsFront,
                                    //~ newFraction
                                //~ );
                            //~ }
                            //~ else
                            //~ {
                                //~ allWeight[nei] = 0.0;
                                //~ allCellTypes[nei] = HOLE;
                                //~ seedCell
                                //~ (
                                    //~ nei,
                                    //~ 1.0,
                                    //~ newIsFront,
                                    //~ newFraction
                                //~ );
                            //~ }
                        //~ }
                    }
                }
            }
        }

        syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
        //~ syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());

        isFront.transfer(newIsFront);
        //~ fraction.transfer(newFraction);
        iPushFront -= 1;
    }

这个实现起来其实很简单，就是用户设置往外挖的层数iPushFront，在每一次循环时，对HOLE外围的cell继续设置成HOLE，由此实现了向外挖洞，而且不影响其他的东西，十分巧妙。

文章版权归作者所有，未经允许请勿转载。

THE END