Just for fun, let’s see what all this actually means in local coordiantes.  It’s a good idea to get our hands dirty too.  To me at least, this makes the operators seem somewhat more friendly.

Choose local coordinates {x^1, \dots, x^n}. Let

\displaystyle g_{ij}(x) = < \partial_i, \partial_j >.

Let {g = \det( g_{ij})}. Note that {g>0}.

Henceforth, we will use the Einstein summation convention: all repeated indices are to be summed over, unless otherwise stated. For instance, the inner product is the 2-tensor

\displaystyle \frac{1}{2} g_{ij} dx^i \otimes dx^j.

Div

Now let {X = X^j \partial_j} be a vector field. We want to compute {\mathrm{div} X} in terms of the quantities {X^i} and {g}. First, it will be necessary to compute the form {dV}. I claim that, if we take {x^i} to be oriented coordinates, then

\displaystyle dV = g^{1/2} dx^1 \wedge \dots \wedge dx^n.

We will prove a more general fact. Whenever {f^i} is an oriented basis for the cotangent space (not necessarily orthonormal), and {g_{ij} = <f^i, f^j>}, {g} is defined as before:

\displaystyle dV = g^{1/2} f^1 \wedge \dots \wedge f^n.

Indeed, this is clear when {f^i} is orthonormal. And it is clear that the above quantity is actually invariant under changes of frame—this is a fact of linear algebra.

So now, with this we can tackle the problem of computing {\mathrm{div} X}. We have:

\displaystyle L_X (dV) = di_X (dV) + i_X d(dV).

I’m not terribly thrilled with this notation. {dV} looks like the exterior differential of some {V}, which it is not. However, this seems to be the norm. Anyhow, {d(dV)} is zero, not because {d^2=0}, but rather because {dV} is an {n}-form. As a result, we are reduced to computing {d i_X(dV)}, which turns out to be rather straightforward. First of all:

\displaystyle i_X(dV) = \sum_j (-1)^{j-1} X^j g^{1/2} dx^1 \wedge \dots \widehat{dx^j} \dots \wedge dx^n,

by definition of how the wedge product works and the interior product. (I guess this is not completely standard. If you use a different definition, you’d get a constant factor added.) Now if we apply {d} to this we find:

\displaystyle \sum_j \partial_j( X^j g^{1/2}) dx^1 \wedge \dots \wedge dx^n ,

and as a result:

\displaystyle \boxed{ \mathrm{div} X = g^{-1/2} \partial_j( X^j g^{1/2}).}

From this the property {\mathrm{div}(Xf) = f \mathrm{div} X + Xf} is clear, since {Xf = X^i \partial_i f}.

Incidentally, to say that {\mathrm{div} X = 0} is to say that the flows of {X} preserve the volume, so it makes sense to call {X} incompressible.

Grad

Now we want to do the same for {\mathrm{grad}}. First, {df = \partial_i f dx^i}. The transformation of a vector {v_i dx^i} to {v^i \partial_i} works by {v^i = g^{ij} v_j}, where {g^{ij}} is the inverse of the matrix {g_{ij}}, I claim. This can be seen as follows:

\displaystyle < g^{ij} v_j \partial_i , w^k \partial_k> = g^{ji} w^k v_j g_{ik} = \delta_{jk} v_j w^k = v_j w^j .

What this rather uninspiring statement means is that { g^{ij} v_j \partial_i} acts as functional on the tangent space via the inner product in the same way that the cotangent vector {v_i dx^i} does. In particular, it proves the claim I just made. (Cf. also the Wikipedia article on raising and lowering indices.) So

\displaystyle \mathrm{grad} f = g^{ij} \partial_j f \partial_i .

The Laplacian

Finally, as a result, we can get the representation of {\Delta f} in local coordinates:

\displaystyle \mathrm{div} \ \mathrm{grad} f = \mathrm{div}( g^{ij} \partial_j f \partial_i ) = g^{-1/2} \partial_i \left( g^{1/2} g^{ij} \partial_j f \right).

In particular, we find that {\Delta} as an operator is {g^{ij} \partial_i \partial_j f} plus something solely first order. Since {g^{ij}} is a positive-definite matrix, this will imply (I will later define this) that the Laplacian is an elliptic operator.