<p>I messed up I am doing a double negative and that is causing the problem:</p>
<p><code>const float laplace_coeff = -0.5f * ((1.f / tan(alpha)) + (1.f / tan(beta)));</code></p>
<p>This should be positive.</p>

CC BY-SA 4.0
2021-04-15T19:45:29.213Z
2021-04-15T19:45:29.213Z
Makogan
<p>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.</p>
<p>The crux of the issue is we are solving a linear system <span class="math-container">$(M - tL_C)u = \delta_y$</span>. In this case <span class="math-container">$M$</span> is a diagonal matrix, <span class="math-container">$L_C$</span> is the discrete Laplace operator, <span class="math-container">$t$</span> is a constant time step equal to some epsilon and <span class="math-container">$\delta_y$</span> is a sparse vector with only one non-zero entry which is equal to 1.</p>
<p>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:</p>
<p><a href="https://i.sstatic.net/j22KI.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/j22KI.png" alt="" /></a></p>
<p>Which looks like random noise to me, as if I messed up the calculation.</p>
<p>I implemented the solver using Eigen, this is the snippet:</p>
<pre><code>Eigen::DiagonalMatrix&lt;float, Eigen::Dynamic&gt; areas(sources.size());
    Eigen::SparseMatrix&lt;float&gt; laplace_operator(sources.size(), sources.size());

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

    areas.setZero();

    for(uint e = 0; e &lt; mesh.Edges().size(); e += 2)
    {
        const auto&amp; 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&amp; v: mesh.Verts())
    {
        std::vector&lt;HMesh&lt;VertexData&gt;::MFace*&gt; faces = v.ContainingFaces();
        for(auto f: faces)
        {
            areas.diagonal()[count] += f-&gt;Area();
        }

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

        count++;
    }

    SparseLU&lt;SparseMatrix&lt;float&gt;, COLAMDOrdering&lt;int&gt; &gt; solver;

    const Eigen::SparseMatrix&lt;float&gt; a =
        Eigen::SparseMatrix&lt;float&gt;(areas) - 0.01 * laplace_operator;
    solver.analyzePattern(a);
    solver.factorize(a);

    return solver.solve(sources);
</code></pre>
<p>I have checked and <code>sources</code> 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.</p>
<p>I suspect my mistake might be creating the laplace operator, but I am not sure.</p>
<p><a href="https://i.sstatic.net/9fusI.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/9fusI.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/WXeUW.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/WXeUW.png" alt="enter image description here" /></a></p>
<p><a href="https://i.sstatic.net/GE9OR.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/GE9OR.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2021-04-15T06:53:47.300Z
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`.
gilgamec
CC BY-SA 4.0
2021-04-15T16:02:25.087Z
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.
Makogan
CC BY-SA 4.0
2021-04-15T16:23:32.557Z
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.
Makogan
<mathematics><computational-geometry><geometry><mesh><linear-algebra>
Computing heat diffusion creates weird results
Makogan
# Computing heat diffusion creates weird results

(Asked by: username_0 on 2021-04-14)

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:

[![](https://i.sstatic.net/j22KI.png)](https://i.sstatic.net/j22KI.png)

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](https://i.sstatic.net/9fusI.png)](https://i.sstatic.net/9fusI.png)

[![enter image description here](https://i.sstatic.net/WXeUW.png)](https://i.sstatic.net/WXeUW.png)

[![enter image description here](https://i.sstatic.net/GE9OR.png)](https://i.sstatic.net/GE9OR.png)

**Tags:** mathematics, computational-geometry, geometry, mesh, linear-algebra

> **Comments:**
> **username_1**: 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`.  
> **username_0**: 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.  
> **username_0**: 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.  

---

## 1 answers

### Answer by username_0 on 2021-04-15. Score: 0 | Accepted answer

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.


---

