Discrete Differential Geometry Lab-1: Discrete Exterior Calculus

Posted on Edited on

每天午饭后做一个小时的几何有助于健胃消食

这次实验的主题是离散曲面上的外微分算子，核心是实现55个离散算子。这是理解曲面外微分运算的关键，也是后续继续学习的基础。

由于我们本次实验讨论的d\text{d}算子和⋆\star算子都是线性算子，并且我们考虑的是它们在离散曲面上的离散版本，维度是有限的，我们可以用矩阵来表达这些算子。实验要求也就是构造出正确的矩阵。

1. Hodge Star of Discrete 0-form\text{1. Hodge Star of Discrete 0-form}

这里要实现离散00形式的霍奇星算子⋆^\hat{\star}。

00形式也就是标量场，其离散版本也就是在离散曲面上给每个点都赋一个权值。由于我们的曲面是一个二维流形，在上面对00形式取霍奇星应该得到22形式，直观理解就是由点到“面”。

在课上，我们得出了离散微分形式的霍奇星算子的一个近似：

⋆^αi^=∣σi⋆∣∣σi∣αi^\hat{\star}\hat{\alpha_i}=\frac{|\sigma_i^\star|}{|\sigma_i|}\hat{\alpha_i}

对于离散微分kk形式α^\hat{\alpha}，如果我们对其施加离散霍奇星算子⋆^\hat{\star}，那么得到的应当是一个离散n−kn-k形式，它们之间系数的数量关系满足上面的等式。

特别的，在00形式的情况下，我们可以认为∣σi∣=1|\sigma_i|=1，而要计算∣σi⋆∣|\sigma_i^\star|，可以使用书上提及的公式——将与点σi\sigma_i相邻的所有面的面积加起来除以33即可。

1
2
3
4
5
6
7
double VertexPositionGeometry::barycentricDualArea(Vertex v) const {
    auto adj = v.adjacentFaces();
    double sum = 0.0;
    for(Face f : adj)
        sum += faceArea(f);
    return sum / 3.0; 
}

接下来构造矩阵就简单了，我们需要的矩阵是一个∣V∣×∣V∣|V|\times |V|的对角矩阵，对角线上的值就是所有的∣σi⋆∣∣σi∣\frac{|\sigma_i^\star|}{|\sigma_i|}。

1
2
3
4
5
6
7
8
9
10
11
12
13
SparseMatrix<double> VertexPositionGeometry::buildHodgeStar0Form() const
{
    Eigen::SparseMatrix<double> spm;
    std::vector<Eigen::Triplet<double>>tList;
    spm.resize(mesh.nVertices(), mesh.nVertices());
    for(Vertex v : mesh.vertices())
    {
        size_t idx = vertexIndices[v];
        tList.emplace_back(idx, idx, barycentricDualArea(v));
    }
    spm.setFromTriplets(tList.begin(), tList.end());
    return spm; // placeholder
}

2. Hodge Star of Discrete 1-form\text{2. Hodge Star of Discrete 1-form}

构造11形式的离散霍奇星矩阵。构造离散霍奇星矩阵的过程其实都大差不差，主要有区别的就是计算∣σi⋆∣∣σi∣\frac{|\sigma_i^\star|}{|\sigma_i|}的过程不一样。

离散11形式对应的是在边上赋权。对于11形式而言，∣σi∣|\sigma_i|就是边长度。而∣σi⋆∣|\sigma_i^\star|是其对偶边的长度，课程中也给出了计算公式——12(cot⁡α+cot⁡β)\frac{1}{2}(\cot\alpha+\cot \beta)，其中α\alpha和β\beta是该边在两个相邻面的对角。

因为我们考虑的是二维流形，所以与一个边相邻的只有两个面，这样这个公式里α\alpha和β\beta就是能唯一确定的。

1
2
3
4
5
6
7
8
9
10
11
double VertexPositionGeometry::cotan(Halfedge he) const {

    Halfedge nhe = he.next();
    Halfedge nnhe = nhe.next();
    assert(nnhe.tipVertex() == he.tailVertex());
    double c = 0.5 * (halfedgeVector(nhe).norm2() + halfedgeVector(nnhe).norm2() - halfedgeVector(he).norm2())
                / (halfedgeVector(nnhe).norm() * halfedgeVector(nhe).norm());
    assert(abs(c) <= 1.0);
    double s = std::sqrt(1.0 - c * c);
    return c / s; // placeholder
}

构造矩阵也非常容易了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
SparseMatrix<double> VertexPositionGeometry::buildHodgeStar1Form() const
{
    Eigen::SparseMatrix<double> spm;
    std::vector<Eigen::Triplet<double>> tList;
    spm.resize(mesh.nEdges(), mesh.nEdges());
    for(Edge e : mesh.edges())
    {
        size_t idx = edgeIndices[e];
        Halfedge he = e.halfedge();
        tList.emplace_back(idx, idx, 0.5 * (cotan(he) + cotan(he.sibling())));
    }
    spm.setFromTriplets(tList.begin(), tList.end());
    return spm; // placeholder
}

3. Hodge Star of Discrete 2-form\text{3. Hodge Star of Discrete 2-form}

这个就更简单了，22形式相当于面，那么∣σi∣|\sigma_i|就是面的面积，并且∣σi⋆∣=1|\sigma_i^\star|=1。

1
2
3
4
5
6
7
8
9
10
11
12
13
SparseMatrix<double> VertexPositionGeometry::buildHodgeStar2Form() const {

    Eigen::SparseMatrix<double> spm;
    std::vector<Eigen::Triplet<double>> tList;
    spm.resize(mesh.nFaces(), mesh.nFaces());
    for(Face f : mesh.faces())
    {
        size_t idx = faceIndices[f];
        tList.emplace_back(idx, idx, 1.0 / faceArea(f));
    }
    spm.setFromTriplets(tList.begin(), tList.end());
    return spm; // placeholder
}

4. Differential Operator of 0-form\text{4. Differential Operator of 0-form}

下面要实现的是0,10,1形式的离散外微分算子d\text{d}。

课上已经讲过如何使用斯托克斯定理推导离散形式的微分形式和外微分算子。对于00形式而言，求外微分得到11形式，相当于由“点权”到“边权”。假设边是ee，那么要求离散的11形式d^α^\hat{d}\hat{\alpha}，稍微算算就知道

d^α^=∫edα=∫∂eα=α(A)−α(B)\hat{\text{d}}\hat{\alpha}=\int_e \text{d}\alpha = \int_{\partial e}\alpha=\alpha(A) - \alpha(B)

对于00形式而言，α^=α\hat{\alpha}=\alpha，所以我们只需要用这条边的两端点的点权作差作为边权即可，写出矩阵也就很容易。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
SparseMatrix<double> VertexPositionGeometry::buildExteriorDerivative0Form() const {
    Eigen::SparseMatrix<double> spm;
    std::vector<Eigen::Triplet<double> >tList;
    spm.resize(mesh.nEdges(), mesh.nVertices());
    for(Edge e : mesh.edges())
    {
        size_t a = vertexIndices[e.firstVertex()];
        size_t b = vertexIndices[e.secondVertex()];
        uint idx = edgeIndices[e];
        tList.emplace_back(idx, a, -1.0);
        tList.emplace_back(idx, b, 1.0);
    }
    spm.setFromTriplets(tList.begin(), tList.end());
    return spm; // placeholder
}

5. Differential Operator of 1-form\text{5. Differential Operator of 1-form}

11形式的外微分得到22形式，也就是由边权计算面权。用类似上面的推导我们可以知道，这个计算也就是将某个面的所有边以某种固定的顺序（逆时针/顺时针）环绕一遍，根据边的定向累加或者减去边的边权。

这里有个问题在于边的定向，“边”本身是无向的，因此我们实际不可能在环绕的时候确定绕行方向和边的方向是否同向。每条边的定向实际上是可以人为规定的，但是仅从面来考虑并不知道边的定向。对于和同一条边相邻的两个面，它们的绕行方向如果是相同的（都为逆时针/顺时针），那么在经过同一条公共边的时候两个面相对于这条边的前进方向必定相反，因此必然有一个面是累加边权，有一个面是减去边权，而我们从面来考虑的话，由于不知道当前的绕行方向（如果我们获取面的所有边那么自然不知道我们是以什么顺序绕行），自然也无法知道应该累加还是减去边权。

所以一个简单的方法是从边来考虑，对于一条边获取其对应的半边，将半边所对应的面的权值增加边权，将另一面减去边权即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
SparseMatrix<double> VertexPositionGeometry::buildExteriorDerivative1Form() const 
{
    Eigen::SparseMatrix<double> spm;
    spm.resize(mesh.nFaces(), mesh.nEdges());
    std::vector<Eigen::Triplet<double>> tList;
    for(Edge e : mesh.edges())
    {
        size_t idx = edgeIndices[e];
        Face fA = e.halfedge().face();
        Face fB = e.halfedge().sibling().face();
        size_t idxA = faceIndices[fA];
        size_t idxB = faceIndices[fB];
        tList.emplace_back(idxA, idx, 1.0);
        tList.emplace_back(idxB, idx, -1.0);
    }
    spm.setFromTriplets(tList.begin(), tList.end());
    return spm; // placeholder
}

另外，在实现上对于本次实验都有个小细节就是怎么获取元素下标。实际上元素都提供了getIndex()方法供你获取下标。但是我一开始眼瞎没注意所以都是用的XXXIndices[]来获取下标，这样也不是不行，但是必须修改主函数，在构造这些算子矩阵之前调用requireXXXIndices()方法。

到此就结束了，总而言之框架相对还是非常完善，数学也很清晰，难度不大。