# Computing heat diffusion creates weird results

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.   • 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. Apr 15, 2021 at 6:53
• 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. Apr 15, 2021 at 16:02
• 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. Apr 15, 2021 at 16:23

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