在OpenFOAM中获取每个网格的Skewness

常用的网格划分软件中关于网格偏度（Skewness）只有统计结果，通常只能看到平均值、最大值及粗略的分布情况，而无法看到每一个网格单元对应的Skewness。因此，这里我希望借助OpenFOAM建立一个生成每个网格单元对应Skewness的小工具。

在OpenFOAM中，可以利用checkMesh工具检查当前网格的质量。checkMesh检测结果包含了网格的Skewness，其输出举例如下：

# ......
Checking geometry...
    Overall domain bounding box (-0.50095 -0.499013 -1) (0.55 0.499013 1.1)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (4.31927e-16 3.08134e-16 1.46623e-16) OK.
    Max cell openness = 2.90386e-16 OK.
    Max aspect ratio = 3.69768 OK.
    Minimum face area = 0.000292504. Maximum face area = 0.00215511.  Face area magnitudes OK.
    Min volume = 8.27995e-06. Max volume = 7.0365e-05.  Total volume = 1.14596.  Cell volumes OK.
    Mesh non-orthogonality Max: 39.9086 average: 4.756
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 3.7248 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

可以看到，checkMesh工具仅仅能够给出最大的Skewness，而无法给出所有网格单元Skewness的详细信息。但无论如何，我们也可以借助checkMesh工具来研究与Skewness计算相关的代码。

checkMesh工具位于$FOAM_APP/utilities/mesh/manipulation/checkMesh文件夹中，该文件夹下包含诸多文件，但通过内容搜索可以知道，与计算Skewness相关的代码仅在checkGeometry.C文件中。相关的内容如下：

faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
if (mesh.checkFaceSkewness(true, &faces))
{
    noFailedChecks++;

    label nFaces = returnReduce(faces.size(), sumOp<label>());

    if (nFaces > 0)
    {
        Info << "  <<Writing " << nFaces
             << " skew faces to set " << faces.name() << endl;
        faces.instance() = mesh.pointsInstance();
        faces.write();
    }
}

这段代码的含义暂且不细研究。联想到checkMesh工具仅能够给出网格整体的最大Skewness，因此我查了一下网格整体Skewness的计算方法。Fluent用户文档中给出：

Skewness is defined as the difference between the shape of the cell and the shape of an equilateral cell of equivalent volume. Highly skewed cells can decrease accuracy and destabilize the solution. For example, optimal quadrilateral meshes will have vertex angles close to 90 degrees, while triangular meshes should preferably have angles of close to 60 degrees and have all angles less than 90 degrees. A general rule is that the maximum skewness for a triangular/tetrahedral mesh in most flows should be kept below 0.95, with an average value that is less than 0.33. A maximum value above 0.95 may lead to convergence difficulties and may require changing the solver controls, such as reducing under-relaxation factors and/or switching to the pressure-based coupled solver.

即Skewness是一个面元与规则面元之间的变形差异情况。此外，Fluent用户文档中也给出了面元的Skewness与网格单元、网格整体的Skewness的计算方式：

Normalized Equiangular Skewness
In the normalized angle deviation method, skewness is defined (in general) as
$$max\left[\frac{{\theta }_{max}-{\theta }_e}{180-{\theta }_e},\frac{{\theta }_e-{\theta }_{min}}{{\theta }_e}\right]$$
where
${\theta }_{max}$ = largest angle in the face or cell
${\theta }_{min}$ = smallest angle in the face or cell
${\theta }_e$ = angle for an equiangular face/cell (such as 60 for a triangle, 90 for a quad, and so on)
The cell skewness will be the maximum skewness computed for any face. For example, an ideal pyramid (skewness=0) is one in which the 4 triangular faces are equilateral (and equiangular) and the quadrilateral base face is a square.

因此，一个网格单元的Skewness值即为它所包含的面元的最大Skewness值。而网格整体的Skewness值即为所有面元的最大Skewness值。显而易见，面元Skewness较大会影响插值精度，因此在绘制网格时，需要注意生成网格的Skewness范围。

接下来我们就寻找checkFaceSkewness()函数所在的位置。由于前面提到，此函数在checkMesh工具中用于对网格的Skewness进行检查，因此其一定包含计算Skewness的过程。在OpenFOAM的源码文件夹$FOAM_SRC/OpenFOAM中，我们通过Ag/Ack工具搜索checkFaceSkewness字串，可以看到包含此函数的文件主要有两个：

• meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
• meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

通过比较这两个源文件中checkFaceSkewness()函数所包含的输出信息我们可以知道，checkMesh工具使用的是primitiveMeshCheck.C中所定义的checkFaceSkewness()。到这里，令人惊喜的事情发生了，在checkFaceSkewness()函数中，可以发现它使用了：

tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness(...)

这个函数来计算某个面的Skewness，因此我们继续寻找这个函数，终于在meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C这个文件中看到了faceSkewness()这个函数的定义：

Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
(
    const primitiveMesh& mesh,
    const pointField& p,
    const vectorField& fCtrs,
    const vectorField& fAreas,
    const vectorField& cellCtrs
)
{
    const labelList& own = mesh.faceOwner();
    const labelList& nei = mesh.faceNeighbour();

    tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
    scalarField& skew = tskew.ref();

    forAll(nei, facei)
    {
        skew[facei] = faceSkewness
        (
            mesh,
            p,
            fCtrs,
            fAreas,

            facei,
            cellCtrs[own[facei]],
            cellCtrs[nei[facei]]
        );
    }


    // Boundary faces: consider them to have only skewness error.
    // (i.e. treat as if mirror cell on other side)

    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
    {
        skew[facei] = boundaryFaceSkewness
        (
            mesh,
            p,
            fCtrs,
            fAreas,
            facei,
            cellCtrs[own[facei]]
        );
    }

    return tskew;
}

该函数在计算一个非壁面面元的Skewness时，将其与包含此面元的两个网格单元一起分析。输入的数据包括该面元对应顶点的坐标、两个网格单元中心连线向量、面元中心与网格单元中心连线向量、面元的面法相向量。而壁面面元也有自己的特别计算方式。通过该函数，可以输出当前面元的Skewness。

根据以上信息，我所写的计算每个网格单元Skewness的小工具步骤包括：

• 读取网格信息
• 使用faceSkewness()函数计算所有面元的Skewness
• 根据面元与网格单元的对应关系，将每个面元的Skewness映射到对应的网格单元上，并只保留最大的Skewness
• 输出网格单元上保存的Skewness

源码如下：

#include "IFstream.H"
#include "OFstream.H"
#include "fvCFD.H"
#include "primitiveMesh.H"
#include "primitiveMeshTools.H"
#include "meshSearch.H"

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    const pointField &p = mesh.points();    // All points
    const vectorField &C = mesh.C();        // All cell centers
    const vectorField &fC = mesh.Cf();      // All face centers
    const vectorField &fA = mesh.Sf();      // All face area vectors

    /* const faceList &faces = mesh.faces(); */
    const labelList &faceOwn = mesh.faceOwner();        // face owner cell
    const labelList &faceNei = mesh.faceNeighbour();    // face neighbour cell

/*  Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
    (
        const primitiveMesh& mesh,
        const pointField& p,
        const vectorField& fCtrs,
        const vectorField& fAreas,
        const vectorField& cellCtrs
    )*/
    // Calculate the skewness of all faces.
    scalarField fskewness =
        primitiveMeshTools::faceSkewness(mesh, p, fC, fA, C);

    // Initialize the skewness of all cells.
    scalarField cskewness = mesh.V();
    forAll(cskewness, cellId) {
        cskewness[cellId] = 0.0;
    }

    // Calculate the skewness of all cells.
    label faceNeiId = 0;
    forAll(fskewness, faceId) {
        cskewness[faceOwn[faceId]]
            = max(cskewness[faceOwn[faceId]], fskewness[faceId]);
        if (mesh.isInternalFace(faceId)) {
            cskewness[faceNei[faceNeiId]]
                = max(cskewness[faceNei[faceNeiId]], fskewness[faceId]);
            faceNeiId++;
        }
    }

    // Print.
    double avg = 0;
    int num = 0;
    forAll(cskewness, cellId) {
        /* Info << cellId << "\t" << cskewness[cellId] << endl; */
        avg += cskewness[cellId];
        num++;
    }

    avg /= num;
    Info << avg << endl;

    return 0;
}