在我之前的博客里最后提到一个解决方案是对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，由此实现了向外挖洞，而且不影响其他的东西，十分巧妙。
暂无评论内容