Imagine a polyhedral, discrete surface embedded in $\mathbb{R}^3$. Its faces are all triangles. For each vertex, one can compute the discrete mean and Gaussian curvatures and evaluate the sum of square principal curvatures, $k_1^2 + k_2^2$.
The goal is to evaluate the bending energy at a vertex $v_0$ by computing this integral over all incident triangles at the $v_0$ vertex:
$$ \int_{S}{k_1^2+k_2^2 dS} $$
The problem is that there's no analytical expression for the $k_i$ curvatures, they're just scalars associated to each $v_k$ vertex of the mesh.
Is there a way to evaluate this particular surface integral? Or, more clearly, given a random triangle in 3D with scalars attached to its vertices, can a notion of surface integral for a function that just interpolates the scalar values over this triangle be established in a consistent way? (then I can just sum the integrals over all triangles and get a numerical estimate).
NOTE
The purpose of computing that surface integral is to assess the importance of a vertex as a salient feature in a local context for a polyhedral surface. The entire concept should be implementable using a computer programming language (hence the more abstract mathematical quirks and impediments may be overlooked up to an acceptable limit).
For a "quick" reference, a good survey of how some continuous differential geometry concepts are adapted in the discrete context, you can consult CalTech's Graphics Group tech report on the issue: http://multires.caltech.edu/pubs/diffGeoOps.pdf
Solution proposition
Let $f(v) = (k_1^2 + k_2^2)(v)$ be a positive scalar function defined on the polyhedral surface. Evaluating $\int_S{f dS}$ can be done by summing all $\int_{\Delta{ijk}}{f dS}$ over all $(v_i,v_j,v_k)$ triangular faces of the surface. Then the problem is reduced to evaluating:
$$ \int_{\Delta{ijk}}{f dS} = \int_s \int_t { f(v(s,t)) \left| \frac{\partial v}{\partial s} \times \frac{\partial v}{\partial t}\right| ds dt },$$ where $v(s,t)$ is a parameterization for the $\Delta_{ijk}$ triangle. For example, $v(s,t) = (1-t)((1-s)v_1 + sv_2) + tv_3$, where $(s,t) \in [0,1]^2$.
Now, we require that $f(v(0,0))=f(v_1)=K_1$, $f(v(1,0))=f(v_2)=K_2$ and $f(v(0,1))=f(v_3)=K_3$, i.e. $f(v_i)=K_i$. Hence $f$ can also be a convex combination of $K_1$, $K_2$ and $K_3$, but this time over a triangular domain that is isometric to $\Delta_{v_1v_2v_3}$.
There should exist a parameterization $v^*(s,t)$ of the $\Delta_{v_1v_2v_3}$ triangle such that $\left| \frac{\partial v}{\partial s} \times \frac{\partial v}{\partial t}\right| = 1$ (is this affirmation valid?).
Since the image of $f(\Delta_{v_1v_2v_3})$ is also a triangle (it is a convex combination of $K_i$s), and since there must be a "canonical" parameterization $v^*(s,t)$ of $\Delta$, then we can compute the resulting double integral as the volume under the $K_1K_2K_3$ triangle over the triangular domain that is isometric to $\Delta_{v_1v_2v_3}$. This solid object is a generalized triangular prism (its base planes are not necessarily parallel anymore).
Is there any obvious flaw in the above solution?