Is there a way to calculate mean curvature of a triangular mesh?
This is an answer based on this Computer Graphics SE indicated in the comments (the one by Nathan Reed).

The math is described in this answer, but in short the curvature is calculated by vertex so:


Get all edges from this vertex
And for each edge compare the projection along the edge of the normals at it extremities
Take the mean of all that


Now, as didn't want to just code a copy of it, I propose a little enhancement:

The calculation is also vertex centered, but we go along each of its surrounding face (need to be triangulated), to calculate a curvature weighted by the angle of this face at the vertex we are calculating around.

Once the sum done over all faces, we mean this sum by the total angle around the vertex.

Making so, the curvature of a large angle will be more important than for a small angle. This is like taking the integral of the curvatures around the vertex.

I think this may be close to the theoretical calculus as it consist of taking the curvature for each cutting plane turning around the vertex normal.

In both algorithms, the base calculus is the same and this corresponds to the CGSE answer:

def curvature_along_edge( vert, other ):
    normal_diff = other.normal - vert.normal
    vert_diff = other.co - vert.co
    return normal_diff.dot( vert_diff ) / vert_diff.length_squared


But in the case we want to use angles, we need to organize a surrounding loop (ring) over the vertex for which we want the curvature, to have the several angles with the good signs.

For instance here:



We want the ring/loop:


For vertex 0: 5, 8, 4
For vertex 8: 1, 6, 3, 7, 2, 4, 0, 5, (1)


etc.

These loops are counterclockwise, as vertices of Blender's polygons are.

This is done by:

# Get vertices in the face order but starting from a given vert
def following_verts_of_vert( vert, face ):
    i0 = index_of( vert, face.verts )
    i1 = (i0 + 1) % 3
    i2 = (i0 + 2) % 3
    return face.verts[i0], face.verts[i1], face.verts[i2]

# Create the oriented ring around vert
def ring_from_vert( vert ):
    vertices = []
    for face in vert.link_faces:
        i0, i1, i2 = following_verts_of_vert( vert, face )
        vertices.append( [i1, i2] )
    result = vertices[0]    
    added = True
    while added and len(vertices):
        added = False
        prev = search_link( result[0], vertices, 1 )
        if prev:
            result = [prev[0]] + result
            vertices.remove( prev )
            added = True
        next = search_link( result[-1], vertices, 0 )
        if next and next[1] not in result:
            result.append( next[1] )
            vertices.remove( next )
            added = True
    return result


(this code above is not really needed, or could be optimized, but was my first though to have a constant orientation reference for the curvature)

So that finally, the mean curvature around a vertex is calculated by:

def angle_between_edges( vert, other1, other2 ):
    edge1 = other1.co - vert.co
    edge2 = other2.co - vert.co
    product = edge1.cross( edge2 )
    sinus = product.length / (edge1.length * edge2.length)
    return asin( min(1.0, sinus) )

def mean_curvature_vert( vert ):
    ring = ring_from_vert( vert )
    ring_curvatures = [curvature_along_edge( vert, other ) for other in ring]
    total_angle = 0.0
    curvature = 0.0
    for i in range(len(ring)-1):
        angle = angle_between_edges( vert, ring[i], ring[i+1] )
        total_angle += angle
        curvature += angle * (ring_curvatures[i] + ring_curvatures[i+1])

    return curvature / (2.0 * total_angle)


Here is the compared results:



They are close, but the comparison I've used needs some explanation: the resulting curvatures are normalized so that they fit to the interval for vertex groups [0, 1]:

In short, this vertex group assignment compares the contrasts but not the values.

This is done so:

def assign_to_vertex_group( obj, group_name, curvatures ):
    vertex_group = ensure_vertex_group( obj, group_name )

    curvatures = [abs(c) for c in curvatures]

    min_curvature = min( curvatures )
    max_curvature = max( curvatures )
    vg_fac = 1.0 / (max_curvature - min_curvature) if max_curvature != min_curvature else 1.0

    for i, vert in enumerate( obj.data.vertices ):
        vertex_group.add( [vert.index], (curvatures[i] - min_curvature) * vg_fac, 'REPLACE' )


Here is the blend file with the 2 scripts named : cgseSimple and cgseSimple2.



Note: 

Also included Mike's implementation (mike text) + the same rewritten a bit (mikesPaper), for understanding and comparison purpose. But from all that, it is hard to know what is the more accurate (how to determinate which is true?). I was not able to understand how this algo is able to determinate curvature orientation (convex vs. concave inclinations).

Also tested all that on big meshes (not included due to the blend file size limit). But that accentuate the diff between the different approaches.