The surface gradient operator is defined as follows
$$\eqalign{
& \mathop \nabla \limits^s = \left( {{\bf{I}} - {\bf{n}} \otimes {\bf{n}}} \right).\nabla \cr
& = {\bf{I}}.\nabla - \left( {{\bf{n}} \otimes {\bf{n}}} \right).\nabla \cr
& = \nabla - \left( {{\bf{n}}.\nabla } \right){\bf{n}} \cr}\tag{1}$$
- $\bf{n}$ is the unit normal vector
- $\bf{I}$ is the second order identity tenor
- $\otimes$ is the tensor product
- $.$ is the scalar product
As you can see in $(1)$ we have subtracted the normal component of the $\nabla $ from it and hence the name surface gradient.
Use $(1)$ to derive your formula. Consider the following
$$\mathop \nabla \limits^s .{\bf{F}} = \left( {\nabla - \left( {{\bf{n}}.\nabla } \right){\bf{n}}} \right).{\bf{F}} = \nabla .{\bf{F}} - \left( {{\bf{n}}.\nabla } \right){\bf{n}}.{\bf{F}} = \nabla .{\bf{F}} - {\bf{n}}.\nabla \left( {{\bf{n}}.{\bf{F}}} \right)\tag{2}$$
Now, if you put ${\bf{F}} = \mathop \nabla \limits^s u$ you can have
$$\mathop \nabla \limits^s .\mathop \nabla \limits^s u = \nabla .\mathop \nabla \limits^s u - {\bf{n}}.\nabla \left( {{\bf{n}}.\mathop \nabla \limits^s u} \right)\tag{3}$$
but
$${\bf{n}}.\mathop \nabla \limits^s u = {\bf{n}}.\left( {\nabla u - \left( {{\bf{n}}.\nabla u} \right){\bf{n}}} \right) = {\bf{n}}.\nabla u - \left( {{\bf{n}}.\nabla u} \right)\left( {{\bf{n}}.{\bf{n}}} \right) = {\bf{n}}.\nabla u - {\bf{n}}.\nabla u = 0\tag{4}$$
and hence
$$\mathop \nabla \limits^s .\mathop \nabla \limits^s u = \nabla .\mathop \nabla \limits^s u\tag{5}$$
I'll try to put something together to help you a bit here:
If you are working with manifolds, you usually only have local descriptions of a manifold $M$, given by coordinate systems $x:M\cap U\rightarrow \mathbb{R}^n$ (here $U$ is some neighbourhood on $M$ diffeomorphis to a ball in Euclidean space).
On a manifold you can do differential calculus by doing it in the image of the coordinate systems. With this approach, differential operators must be shown to be well defined, which basicall means that they behave in the same way regardless of the choice of coordinate system.
If you, in addition, have a Riemannian metric on the manifolds, several constructions known from differential calculus carry over to manifolds. One of these is the Laplacian, which is defined as the trace of the Hessian (read: second derivative, whatever this means precisely). It turns out that the resulting expression of the Laplacian (also called Laplace Beltrami Operator) is given by the expression you have written down -- there is a typo in there, the formula shoud read
$$\Delta f = \frac{1}{\sqrt{G}}\sum_{i,j} \frac{\partial}{\partial x^i}\left( \sqrt{G}g^{ij}\frac{\partial f}{\partial x^j }\right)$$
Here $g^{ij}$ is the inverse of the metric tensor $g_{ij} =\langle \frac{\partial}{\partial x^i}, \frac{\partial}{\partial x^j}\rangle $
As mentioned earlier, usually you have this represention only locally, and it looks rather complicated. In some cases you can find an easier reprentation (e.g. when the metric tensor is diagonal, or even better, a function times the identity). In the one dimensional case you have a curve and the metric tensor is just a function. If you choose the inverse of an arclength parametrization as coordinate system, the metric tensor becomes (as you indicated in your post) a constant, and the operator reduces to the ordinary second derivative. In case of the circle you can choose a parametrization which covers all of the manifold with the exception of a point, which allows you to identify functions on the circle with periodic functions on the interval. If you put all this together the equation $\Delta f = \lambda f$ reduces to $f^{\prime\prime} = \lambda f$. (The sign of the operator varies among authors).
That's about it. In order to grasp the details and fully understand them in a general setup you will have to dig into Differential Geometry.
Best Answer
Given a point $p\in M$, you can recover the metric at $p$ as follows. Choose any local coordinates $(x^1,\dots,x^n)$ such that $p$ has the coordinate representation $(0,\dots,0)$. For any indices $k,l$, let $f_{kl}(x) = \tfrac 1 2 x^k x^l$. If you expand out the expression for $\Delta (f_{kl})$ and use the fact that $x^k=x^l=0$ at $p$, you'll find that $\Delta(f_{kl})(p) = g^{kl}$. Then you can use matrix inversion to recover $g_{kl}$.
If your oracle is restricted to acting only on globally defined functions, you can use a bump function to extend the coordinate functions to global smooth functions $(u^1,\dots,u^n)$ that agree with $(x^1,\dots,x^n)$ in a neighborhood of $p$, and let $f_{kl}(x) = u^k u^l$.
This is a special case of a much more general construction: if $P$ is an $m$th order scalar linear partial differential operator, the coefficients of its highest-order terms determine a coordinate-independent symmetric contravariant $m$-tensor field called the principal symbol, which can be evaluated at any point by applying $P$ to $m$-fold products of coordinate functions. The principal symbol of $\Delta$ is the associated cometric, which is the induced metric on $T^*M$.