Why don't we see the covariant derivative in classical mechanics? I am wondering why I have seen the covariant derivative for the first time in general relativity.
Starting from the point that the covariant derivative generalise the concept of derivative in curved space (even if think it is better to consider it as the extension of the derivative such that it is covariant under a change of coordinates). To do that we introduce the Christoffel symbols $\Gamma^i_{jk}$.
In curved space-time we have globally non vanishing Christoffel symbols $\Gamma^i_{jk} \neq 0$, but in general $\Gamma^i_{jk} \neq 0$ doesn't mean we are in curved space-time. For example, if I consider Minkowski space-time with Cartesian coordinates than, thanks to the Lorentz transformation, if the Gammas are zero in a reference frame they are zero in every reference of frame, but I could have $\Gamma^i_{jk} \neq 0$ even in flat space time with polar coordinates, as the Gammas don't transform like a tensor in this case due to the non tensorial part of the transformation law for $\Gamma^i_{jk}$ under a change of basis.
If what I said before is true (a big if), then I'd interpret this in classical mechanics saying that in Cartesian coordinates, the basis vectors {$\hat{e}_x,\hat{e}_y$}, solid to a point of a curve, are constant if the point is moved along the curve.
While I think I can't say the same for {$\hat{e}_r, \hat{e}_{\theta}$}, as moving a point along the curve in this case the tangent vectors to the coordinate lines aren't constant (they rotate while the point is moving). This is why I think I should see the Christoffel symbols even in classical mechanics, to reflect the property of the vectors {$\hat{e}_r, \hat{e}_{\theta}$} that vary along the curve.
 A: The Chrisoffel symbol - or the connection to the metric - or just the connection - is the result of taking the derivative of a vector field - which may cause the resulting vector field to rotate.
To determine if a manifold is intrinsically or extrinsically curved, you need to calculate the Riemann curvature tensor. 
For instance, for Euclidean and Minkowski spaces, the Riemann curvature tensor is zero since both of those spaces are extrinsically flat - or just flat spaces.
However, it is possible to embed an intrinsically curved surface in a flat space - in which case one or possibly more Chrisoffel symbols may not be zero - but the Riemann tensor will still be zero. 
The magic of semi-Riemann manifolds is the connection known as the Levi-Civita which is unique.
A: You don't see the covariant derivative as often because flat space has isometries that make Cartesian coordinates better, and in these coordinates there are no Christoffel symbols, so we use them as much as possible. But look at the formula for the divergence of a function $\mathbf{F} = F^\hat{r} \hat{\mathbf{r}} + F^\hat{\theta} \hat{\mathbf{\theta}}$ in polar coordinates:
$$\nabla \cdot \mathbf{F} = \frac{1}{r} \frac{\partial(r F^\hat{r})}{\partial r} + \frac{1}{r} \frac{\partial F^\hat{\theta}}{\partial \theta} = \frac{\partial F^\hat{r}}{\partial r} + \frac{1}{r} F^\hat{r} + \frac{1}{r} \frac{\partial F^\hat{\theta}}{\partial \theta}.$$
That $1/r$ in the middle term with no derivatives comes from the Christoffel symbols! So the covariant derivative is definitely there, but instead of using the Christoffel symbols, we usually calculate it using the chain rule and the fact that the cartesian basis vectors have zero derivative. The derivatives of the basis vector are after all the Christoffel symbols, so the method is not that different.
One final comment: the orthonormal basis vectors $\{\hat{\mathbf{r}}, \hat{\theta}\}$ in polar coordinates are not the basis vectors $\{\partial/\partial r, \partial/\partial\theta\}$ we know from differential geometry, because the latter are not orthonormal. The relation is simple:
$$\begin{aligned}
\hat{\mathbf{r}} &= \frac{\partial}{\partial r} \\
\hat{\theta} &= \frac{1}{r} \frac{\partial}{\partial\theta},
\end{aligned}$$
so keep this in mind when applying the formulas. In differential geometry we tend to write the components of vectors with respect to the derivative basis, but the formulas we know from more basic calculus (like my divergence formula) are written in terms of the orthonormal basis.
A: Another point to consider is that in Hamiltonian mechanics the symplectic structure is independent of a metric. In the regular, non-degenerate case, this structure may be pulled-back to the tangent bundle and the domain of the Lagrangian formulation. 
Therefore, you need not start at the covariant derivative for classical mechanics and instead may recover a more general, abstract description. 
