Area gradient, mean curvature normal, and the "cotan formula" 
Show that when we take the gradient of the total surface area with respect to the position of one of the vertices, we get
$$\nabla_{f_{i}}\sum_{ijk\in F}A_{ijk}=\frac{1}{2}\sum_{ij\in E}\left(\cot\alpha_{ij}+\cot\beta_{ij}\right)\left(f_{i}-f_{j}\right)$$
where $\alpha_{ij}$ and $\beta_{ij}$ are the two interior angles opposite edge $ij$. The right-hand-side is a discretization of the mean curvature normal $HN$ at vertex $v_{i}$.

(from assignment 2 of Keenan Crane's spring 2021 discrete differential geometry course)
I've been poking at this for a while and haven't found a decent direction to go to get there.  The closest I've gotten is this:

if we have $h_1$, we can use that and $\alpha$ to get the length of the base of the outer portions of each of the two triangles $\ell_1 = h_1 \cot \alpha$; the area is then $\frac{1}{2} \ell_1|h_1|$, and the area gradient is then $\frac{1}{2}\ell_1\hat h_1=\frac{1}{2}h_1\cot \alpha$ (because we don't actually move $\ell_1$ when moving the common vertex), and then the other triangle works the same way, so we get $\frac{1}{2}(h_1\cot \alpha + h_2\cot\beta)$.
However, this doesn't seem to get me anything corresponding to the edge in between - in particular because this isn't truly flat: as long as there's at least some dihedral angle between the two triangles, the two $h$s and the edge won't be coplanar.
What am I missing here?
 A: I figured it out.
The trick is to present the base of the altitude as a linear interpolation between the two base vertices.

We have the triangle's area as $A=\frac{1}{2}(\ell_1 +\ell_2)||C-P||$ and area gradient $\nabla A_C=\frac{1}{2}(\ell_1+\ell_2)\frac{C-P}{||C-P||}$.  Since $D$ is on $AB$ we can write it as $$D=\frac{A\ell_2+B\ell_1}{\ell_1+\ell_2}$$.  Then we can substitute this back into the top of $\nabla A$ to get:
$$
\begin{alignedat}{1}\nabla A_{C} & =\frac{1}{2}(\ell_{1}+\ell_{2})\frac{C-\frac{A\ell_{2}+B\ell_{1}}{\ell_{1}+\ell_{2}}}{||C-D||}\\
 & =\frac{1}{2}(\ell_{1}+\ell_{2})\frac{\frac{C(\ell_{1}+\ell_{2})}{\ell_{1}+\ell_{2}}-\frac{A\ell_{2}+B\ell_{1}}{\ell_{1}+\ell_{2}}}{||C-D||}\\
 & =\frac{1}{2}(\ell_{1}+\ell_{2})\frac{C\ell_{1}+C\ell_{2}-A\ell_{2}-B\ell_{1}}{(\ell_{1}+\ell_{2})||C-D||}\\
 & =\frac{1}{2}\frac{C\ell_{1}+C\ell_{2}-A\ell_{2}-B\ell_{1}}{||C-D||}\\
 & =\frac{1}{2}\left(\frac{(C-A)\ell_{2}}{||C-D||}+\frac{(C-B)\ell_{1}}{||C-D||}\right)\\
 & =\frac{1}{2}\left((C-A)\cot\beta_{CA}+(C-B)\cot\alpha_{CB}\right)
\end{alignedat}$$
And then we can add up all the triangles around the moving vertex and factor the result by taking parts with the same edge in question to get the original result: each edge contributes its vector towards the moving point times the cotangents of its opposing angles to the area gradient.
