numerically integrating a function over triangular elements

I had to do this problem for some C++ code I'm writing for work and so decided to put the result here.

I have a scalar function defined over an irregular 2-dimensional mesh of triangular elements and I need to integrate the function over the surface, for example to find the average value. All I have are the coordinates of the elements.

triangular mesh example, from University of Vermont

Consider an element with points p₀, p₁, and p₂. I can define side vectors r₀ ≡ p₁ - p₀, and r₁ = p₂ - p₁. Then I can determine the area of the triangle a = | r₀ × r₁|.

The integral of the function over the area is the multiplication of the average value of the function times the area. I have the area so I need the average value.

This is of course ill-defined because I don't know how the function varies between the three points. But I can guess that it varies linearly. A linear function in two dimensions is defined by three degrees of freedom, for example an intercept and two slopes, one for each dimension. So I need three constraints to determine a linearization of the function in two dimensions. As long as the three points are not colinear (if they are it doesn't matter because the area of my element is zero anyway) then the three points fully determine this linearization. If it was a 4-sided element there'd be a bit of a problem because then I'd be potentially overconstraining the linearization. But for triangular elements with three points the problem is well defined.

So I should be able to solve this analytically.

Consider two cases: one is where the triangle is equilateral. In this case each point is equidistant from the other two and the average is trivially the average of the three values. Nothing differentiates one point from the others so they must be counted equally.

Here's another case: suppose two points are much closer together than the other, "virtually touching". In this case the average value is an average of the third point (the one further away) and the additional average of the two near points. So it is the weighted sum of 25% of the close points and 50% of the further point.

I can capture these two special cases as follows:

1. Calculate the centroid of the triangle. This is the average of the three positions.
2. Weight the value of each point proportional to the distance from this centroid

This is a bit counter-intuitive and indeed I got it wrong at first because I assumed that the points closer to the center would have more influence over the average value. But this is backwards. The center is drawn to where points are more closely clumped and so I need to de-weight the points in the clump to avoid that portion of the scalar field not be over-represented.

So the result for the average of the function f, with values at points p₀, p₁, and p₂ of f₀, p₁, and p₂, is:

<f_avg> =
   [ f₀ |p₀ - p_avg| + f₁ |p₁ - p_avg| + f₂ |p₂ - p_avg| ] /
   [ |p₀ - p_avg| + |p₁ - p_avg| + |p₂ - p_avg| ]

Is this the most elegant representation Almost certainly it isn't. But it seems to me it's correct, and computers are good at calculating this sort of thing, so I'll let them do their work.

For a 4-sided element, I can divide it into triangles (for example, points 0, 1, and 2, then points 0, 2, and 3), and take the integral of each and sum. Note I may need to be careful here with how I treat my areas, because sub-triangles may have negative area, and to treat them with an absolute value will overestimate the net integral or net area.