On integral of a function over a simplex Help w/the following general calculation and references would be appreciated.
Let $ABC$ be a triangle in the plane. 
 Then for any linear function of two variables $u$. 
$$
\int_{\triangle}|\nabla u|^2=\gamma_{AB}(u(A)-u(B))^2+
\gamma_{AC}(u(A)-u(C))^2+\gamma_{BC}(u(B)-u(C))^2,
$$
 where
$$
\gamma_{AB}=\frac{1}{2}\cot(\angle C),
\gamma_{AC}=\frac{1}{2}\cot(\angle B),
\gamma_{BC}=\frac{1}{2}\cot(\angle A).
$$
What is a good reference for the formula? Is it due to R. Duffin?
Is there generalization to linear functions of three variables? The number of parameters fits (in any $nD$), and it seems like one needs to calculate 6x6 Cayley-Menger like determinants, but it's difficult and the geometric interpretation is not clear.
Thank you.
 A: The solution by @Shuhao Cao is nice. However it is too long to follow. There is a short way to do this. It is easy to check that $\int_{\Delta}|\nabla u|^2$ is invariant under translation  and rotation 
$$ x\to x+h,y\to y+k, x\to x\cos\theta-y\sin\theta,y\to x\sin\theta+y\cos\theta.$$
Thus we can take the coordinates of $A,B,C$ as
$$ A(0,0), B(x_B,0), C(x_C,y_C) $$
respectively, from which we obtain
$$ \cot A=\frac{x_C}{y_C}, \cot B=\frac{x_B-x_C}{y_C}, |\Delta|=\frac{1}{2}x_By_C. $$
From the identity 
$$ \cot A\cot B+\cot B\cot C+\cot C\cot A=1$$
we have
$$ \cot C= \frac{-x_Bx_C+x_C^2+y_C^2}{x_By_C}$$
Let $u=\alpha x+\beta y+\gamma$ and hence $|\nabla u|^2=\alpha^2+\beta^2$. So we have to show the following identity
$$ (\alpha^2+\beta^2)|\Delta|=\frac{1}{2}\alpha^2 x_B^2\cot C+\frac{1}{2}(\alpha x_C-\beta y_C)^2\cot B+\frac{1}{2}[\alpha (x_B-x_C)-\beta y_C]^2\cot A.$$
Note 
\begin{eqnarray*}
RHS&=&\frac{\alpha^2}{2}[x_B^2\cot C+x_C^2\cot B+(x_B-x_C)^2\cot A]\\
&&+\frac{\alpha^2}{2}y_C^2(\cot B+\cot A)+\alpha\beta[-x_Cy_C\cot B+(x_B-x_C)\cot A].
\end{eqnarray*}
It is easy to check
$$ y_C^2(\cot B+\cot C)=x_By_C,-x_Cy_C\cot B+(x_B-x_C)\cot A=0 $$
and 
$$ x_B^2\cot C+x_C^2\cot B+(x_B-x_C)^2\cot A=x_By_C. $$
Thus 
$$ RHS=(\alpha^2+\beta^2)|\Delta|=LHS. $$
A: Let $(x_k, y_k,z_k)$ be four points in 3D. Then for a linear function 
$$u(x,y,z)=Ax+By+Cz+D,$$ 
and
$$ x_{kl}=x_k-x_l, \dots $$
matching coefficients in the identity in question, one gets the following system of equations:
$$
\begin{pmatrix}
x_{12}^2 & x_{13}^2 & \dots & x_{34}^2 \\
y_{12}^2 & y_{13}^2 & \dots & y_{34}^2 \\
z_{12}^2 & \dots & \dots & \dots \\
x_{12}y_{12} & x_{13}y_{13} & \dots & \dots \\
x_{12}z_{12} & \dots & \dots & \dots \\
y_{12}z_{12} & y_{13}z_{13} & \dots & y_{34}z_{34} \\
\end{pmatrix}
\begin{pmatrix}
\gamma_{12} \\ \gamma_{13} \\ \dots \\ \dots \\ \dots \\ \gamma_{34}
\end{pmatrix} =
Vol_{\triangle} 
\begin{pmatrix}
1 \\ 1 \\ 1  \\ 0 \\ 0 \\ 0
\end{pmatrix},
$$
where
$$ x_{kl}+x_{lm}+x_{mk}=0, \dots $$
 Can these expression be simplified? Is there geometric interpretation of the solution?
