I tried to solve the 2D heat equation with the finite element method, using triangles as elements. Currently generated by a Delaunay triangulation. The base function I'm currently using is basically just the first barycentric coordinate. Assuming a triangle with $a, b, c$ as vertices: $$ f(x, y) = \frac{(b_y - c_y) (x - c_x) + (c_x - b_x) (y - c_y)}{(b_y - c_y) (a_x - c_x) + (c_x - b_x) (a_y - c_y)} $$ For integrating the equation I'm using a Gauss4x4 from https://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/quadrature_rules_tri.html over the unit triangle. To transform the triangle($u, v$ are coordinates on the unit triangle): $$ p = a + (b - a) u + (c - a) v $$ and its determinant: $$ (b_x - a_x) (c_y - a_y) - (c_x - a_x) (b_y - a_y) $$ And the weak form of the equation, if the integral over the boundary is zero: $$ -\int_{\Omega} (\frac{\partial u}{\partial x}+\frac{\partial u}{\partial y}) (\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y})= \int_{\Omega} \frac{\partial u}{\partial t}v $$ Now the problem I'm having is that starting from a "heat" patch in the middle the heat flows to the upper left and lower right corner first. Below is picture from the initial state and later when the heat flows into the corners, forming in the long run a diagonal band.

initial heat distribution

enter image description here

I expected the heat to spread in all directions equally. Looking more uniform and like a circle.

Without me posting too much code, may there be an explanation why this is happening?

  • $\begingroup$ The weak form you have written is incorrect. The LHS should be $-\int_\Omega \nabla u\cdot \nabla v \ dx$. If this isn't a typo, this may be causing the issue because the weak form you wrote isn't rotationally invariant and will diffuse in different directions unequally. I'm not putting this as an answer because I'm not sure if it's a typo or not. I can provide more details if you confirm this. $\endgroup$
    – whpowell96
    May 20, 2023 at 22:15
  • $\begingroup$ Correct me if I'm wrong, but that weak form seems to translate to $\partial^2_1 u + 2\partial_{12}^2 u + \partial_2^2 u = \partial_t u$ rather than $\Delta u = \partial_t u$. In cleaner terms, $\nabla \cdot \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} \nabla u$ rather than $\nabla \cdot \nabla u$. $\endgroup$
    – Sardine
    May 20, 2023 at 22:22
  • $\begingroup$ What would be the correct definition of the LHS then, in terms of partial derivatives? I'm pretty sure that this not a typo but rather my mistake $\endgroup$
    – Boiler4562
    May 20, 2023 at 22:28
  • $\begingroup$ $-\int u_xv_x + u_yv_y \ dxdy$ $\endgroup$
    – whpowell96
    May 20, 2023 at 22:30
  • 1
    $\begingroup$ You should write a short answer to this question so it can be resolved. $\endgroup$
    – whpowell96
    May 20, 2023 at 22:48

1 Answer 1


The problem was that the LHS of the weak form was wrong, the correct one is: $$ -\int u_{x} v_{x} + u_{y} v_{y} dxdy $$ Instead of $$ -\int (u_{x}+u_{y}) (v_{x}+v_{y})dxdy $$ Thanks to whpowell96 for spotting the mistake


Your Answer

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

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