@mohui 这个isoAdvector我大概看过，每个有自由液面的网格(0<alpha<1)的等值面不一定一样。根据每个网格中心的alpha值先插值得到网格节点的alpha值，再利用节点的值构建等值面，使等值面以下的水的几何体积等于该网格中alpha*V的值。得到等值面的信息，包括面心位置，法相和速度，可以得到在dt时间内，等值面以下网格面浸水面积随时间的变化，再由已知的phi值计算dt时间内水的通量。也就是在计算这个通量的时候，考虑了等值面的运动。为了确保这个算法的守恒还用了bounding，大概就是如果某个网格接收水以后，其alpha值大于1，多余的部分不是直接砍掉，而是转移的别的网格。

@Tens @李东岳 感谢各位的回复，我们的snap选择了true，也尝试过给SurfaceFeature加密，但是都没有效果，如下图所示。

通过分析，感觉是因为红色部分网格和灰色部分网格处于同一平面内，SHM就识别不了边界。

我们尝试了让红色圆形的部分凸出来一点，或者凹下去一点，但是需要在灰色的网格部分再画一个过度区，这样圆形的边界就能准确的描述出来，但是这样就破坏了原有的几何形状。

各位同仁，最近我们在用snappy画网格的时候遇到一个问题，跟大家请教一下。我们想在平板网格上加几个圆形的压力测点，如下图所示。这是用犀牛建的几何模型。

我们想在用snappy生成网格，并且把黄色圆孔设为独立的边界条件，方便压力的测量。但是生成的网格如下图所示，圆形边界的边缘不能被准确的描述出来，而是呈台阶状。

试了好多种方法，包括边界处加密和调整snappy的参数等等，都不好用，还请大家帮忙看看有没有什么好的办法。多谢！

https://openfoamwiki.net/index.php/Installation/Linux/OpenFOAM-dev/Ubuntu 可以参考第九条 (Save an alias in the personal .bashrc file, simply by running the following command)

@mohui 不是大神，不是大神，互相交流:laughing: 。

1. FTM没了解过，这个isoAdevtor里面，我的理解是，比如二维的正交网格，每个网格节点都有四个网格相连，这个网格节点的值是由周围四个网格中心的值差值而来。
2. 是以alphaV为条件，找到对应的等值面。这个构建等值面的过程是根据每个网格节点处的alpha值，比如二维问题，我们就有四个alpha值，这样可以构首先建四个等值面。这四个等值面对应着四个alphaV，判断一下哪两个alphaV中间包含着我们的目标alphaV，然后就在这个范围内构建一个等值面以下的体积关于等值面的值的函数V=func(isovalue)，这样就可以找到目标等值面。这里面都是纯几何操作。

@mohui 这个isoAdvector我大概看过，每个有自由液面的网格(0<alpha<1)的等值面不一定一样。根据每个网格中心的alpha值先插值得到网格节点的alpha值，再利用节点的值构建等值面，使等值面以下的水的几何体积等于该网格中alpha*V的值。得到等值面的信息，包括面心位置，法相和速度，可以得到在dt时间内，等值面以下网格面浸水面积随时间的变化，再由已知的phi值计算dt时间内水的通量。也就是在计算这个通量的时候，考虑了等值面的运动。为了确保这个算法的守恒还用了bounding，大概就是如果某个网格接收水以后，其alpha值大于1，多余的部分不是直接砍掉，而是转移的别的网格。

@wwzhao @金石为开
我这两天看了这篇文章，里面写的interfoam的理论部分还比较详细。

Evaluating the performance of the two-phase flow solver interFoam

http://iopscience.iop.org/article/10.1088/1749-4699/5/1/014016/meta

@vivian @cfd-china 看来我来晚了:happy:

@wwzhao 谢谢您的回复！

我找到了这个帖子

http://www.cfd-online.com/Forums/openfoam-programming-development/72456-matrix-storage-diagonal-coeffients-calculation.html

但是也没说清楚。不过系数矩阵是对称的话，确实不影响结果。

为了进一步测试，建立一个3X3的网格，并创建一个fvMatrix。

        fvVectorMatrix UEqn_lap
        (
          fvm::laplacian(nu, U)
        );

网格编号如下：

|6|7|8|
|3|4|5|
|0|1|2|

输出 upper 和lower 的address，没有问题

        Info<< "UEqn_lap upper add before sum   " << UEqn_lap.lduAddr().upperAddr() << nl << endl;
        Info<< "UEqn_lap lower add before sum   " << UEqn_lap.lduAddr().lowerAddr() << nl << endl;
UEqn_lap upper add before sum
12
(
1
3
2
4
5
4
6
5
7
8
7
8
)


UEqn_lap lower add before sum
12
(
0
0
1
1
2
3
3
4
4
5
6
7
)

下面把upper lower 和diag重新赋值 并输出

        UEqn_lap.upper() = 100;
        UEqn_lap.diag() = 0;
        UEqn_lap.lower() = 1;

        Info<< "UEqn_lap diag after change   " << UEqn_lap.diag() << nl << endl;
        Info<< "UEqn_lap upper after change   " << UEqn_lap.upper() << nl << endl;
        Info<< "UEqn_lap lower after change   " << UEqn_lap.lower() << nl << endl;

输出为：

UEqn_lap diag after change   9{0}

UEqn_lap upper after change   12{100}

UEqn_lap lower after change   12{1}

调用

        UEqn_lap.negSumDiag();

        Info<< "UEqn_lap diag after negSumDiag   " << UEqn_lap.diag() << nl << endl;

输出结果为：

UEqn_lap diag after negSumDiag   9(-2 -102 -101 -102 -202 -201 -101 -201 -200)

我特意把upper里面的数值设的很大，从输出的结果看，diag 是把lower里面的数加和了，也就是把列中的数加和了。

这个结果让我有点儿弄不明白。

在看关于negSumDiag()的代码的时候遇到一个问题。diag上面的项是由该“行”上面的的off-diag的项相加后取负值得到。为什么下面的代码感觉像是把该“列”的数相加了呢？

template<class Type, class DType, class LUType>
void Foam::LduMatrix<Type, DType, LUType>::negSumDiag()
{
    const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
    const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
    Field<DType>& Diag = diag();

    const unallocLabelList& l = lduAddr().lowerAddr();
    const unallocLabelList& u = lduAddr().upperAddr();

    for (label face=0; face<l.size(); face++)
    {
        Diag[l[face]] -= Lower[face];
        Diag[u[face]] -= Upper[face];
    }
}
 

以这个网格为例

lowerAdd 为 (0 1 2 1 2 3 4 3 )
upperAdd 为 (2 2 3 6 4 5 6 5 )

矩阵的形式为

那么Diag[0] 应该是由 Upper[1]得到。这个和上面的代码不符。如果矩阵是对称的，这个计算没有问题。而如果是 非对称的，是不是就有问题了？

谢谢大家！