Express your problem differently:
If you keep in mind that "the gradient of a function at a point is perpendicular to the level set of the function at that point" - ie that for a function defining a surface, the gradient will be normal to this surface (Q: Why is the gradient normal?) - then all you want to do with your mesh is to compute its vertices' unit normals.
Using this wording, you should find plenty of methods to do so. For instance, here is a simple algorithm described by Perlin here:
Initialize all the vertex normals to the null vector. In the next step, you're going to gradually accumulate values into these vectors.
Loop through all the faces of your mesh. For each face, consider all pairs of successive edges around the face. For example, if your
face has vertices A,B,C,D, in that order, then you would consider the
four pairs of successive edges:
- A→B and B→C
- B→C and C→D
- C→D and D→A
- D→A and A→B
For each pair of successive edges $e_1$ and $e_2$, take the cross product $e_1 × e_2$. Add the resulting vector to the normal of every
vertex that adjoins the face (for example, in this case you would add
to the normals of vertices A,B,C,D).
Loop through all the vertices, and normalize the length of the normal vector of each one.