0
$\begingroup$

I am trying to model heat diffusion on the surface of a mesh. I annexed the most important bits of theory about this topic as screenshots on the question see the bottom.

The crux of the issue is we are solving a linear system $(M - tL_C)u = \delta_y$. In this case $M$ is a diagonal matrix, $L_C$ is the discrete Laplace operator, $t$ is a constant time step equal to some epsilon and $\delta_y$ is a sparse vector with only one non-zero entry which is equal to 1.

My expectation is that, upon solving this system and drawing the resulting vertex distances, I get a fading from white to black starting at the source vertex and slowly fading away as the vertices go away from the source. I am however seeing this:

Which looks like random noise to me, as if I messed up the calculation.

I implemented the solver using Eigen, this is the snippet:

Eigen::DiagonalMatrix<float, Eigen::Dynamic> areas(sources.size());
    Eigen::SparseMatrix<float> laplace_operator(sources.size(), sources.size());

    for(uint i = 0; i < sources.size(); i++)
        laplace_operator.insert(i, i) = 0;

    areas.setZero();

    for(uint e = 0; e < mesh.Edges().size(); e += 2)
    {
        const auto& edge = mesh.Edges()[e];
        const uint i = edge.Vert().ID();
        const uint j = edge.Pair().Vert().ID();

        const Eigen::Vector3f e1 = -edge.Next().Dir().normalized();
        const Eigen::Vector3f e2 = edge.Prev().Dir().normalized();

        const float alpha = acos(e1.dot(e2));

        const Eigen::Vector3f p1 = -edge.Pair().Next().Dir().normalized();
        const Eigen::Vector3f p2 = edge.Pair().Prev().Dir().normalized();

        const float beta = acos(p1.dot(p2));

        const float laplace_coeff = -0.5f * ((1.f / tan(alpha)) + (1.f / tan(beta)));

        laplace_operator.insert(i, j) = laplace_coeff;
        laplace_operator.insert(j, i) = laplace_coeff;

        laplace_operator.coeffRef(i, i) -= laplace_coeff;
        laplace_operator.coeffRef(j, j) -= laplace_coeff;
    }

    uint count = 0;
    for(auto& v: mesh.Verts())
    {
        std::vector<HMesh<VertexData>::MFace*> faces = v.ContainingFaces();
        for(auto f: faces)
        {
            areas.diagonal()[count] += f->Area();
        }

        areas.diagonal()[count] /= 3.f;

        count++;
    }

    SparseLU<SparseMatrix<float>, COLAMDOrdering<int> > solver;

    const Eigen::SparseMatrix<float> a =
        Eigen::SparseMatrix<float>(areas) - 0.01 * laplace_operator;
    solver.analyzePattern(a);
    solver.factorize(a);

    return solver.solve(sources);

I have checked and sources is indeed a sparse vector with a singular source vertex set to 1 so the input is correct. The mesh operations that fetch values are very tested and they are returning what is expected, so I know the error is in the logic of this function.

I suspect my mistake might be creating the laplace operator, but I am not sure.

enter image description here

enter image description here

enter image description here

$\endgroup$
3
  • $\begingroup$ I think you may be computing the wrong cotans. From your figure, if edge is from i to j and edges are ordered clockwise, then $\alpha$ is between edge.Next() and edge.Pair().Prev() -- in your code, that'd be e1 and p2. Similarly, $\beta$ is between e2 and p1. $\endgroup$
    – gilgamec
    Apr 15, 2021 at 6:53
  • $\begingroup$ Are you sure? we are talking in terms of a half edge data structure. edge.Pair().Prev() and edge.Next() should not be on the same face, so i fail to see how they are the correct angles, but perhaps I am just confused here. $\endgroup$
    – Makogan
    Apr 15, 2021 at 16:02
  • $\begingroup$ I am a bit confused. The mesh is a half edge representation. thus edge and its pair correspond to the 2 different faces in the diagram. Basically edge.Next() and edge.Pair().Prev() should correspond to different faces entirely. But maybe I am not seeing something. $\endgroup$
    – Makogan
    Apr 15, 2021 at 16:23

1 Answer 1

0
$\begingroup$

I messed up I am doing a double negative and that is causing the problem:

const float laplace_coeff = -0.5f * ((1.f / tan(alpha)) + (1.f / tan(beta)));

This should be positive.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.