As suggested in @tmonteil 's answer, all necessary material is in the manifolds part, but it needs to be made explicit for the end user. For instance, there should be a function `EuclideanSpace`

so that

```
sage: E = EuclideanSpace(3)
```

creates the 3-dimensional Euclidean space as a smooth manifold endowed with a Riemannian flat metric. Then one could do things like

```
sage: v = E.vector_field([-x*z, x+z, y^2])
sage: c = v.curl()
sage: d = v.div()
sage: spher = E.spherical_frame() # orthonormal frame associated with spherical coordinates
sage: c.display(spher)
...
```

This was on my todo list for a while... @quantum_leopard 's question acts as an efficient reminder: I'm on it ;-)
I will report here (with a Trac ticket number) when the code is ready.

EDIT (18 March 2018): the code implementing vector calculus on Euclidean space is ready for review at Trac #24623.

EDIT (10 May 2018): the ticket #24623 has received a positive review and has been merged in SageMath 8.3.beta0. Alas, it was too late to be integrated in SageMath 8.2, which has been released 5 days ago. However, among the new features of SageMath 8.2, there are the grad, div, curl, etc. operators on a generic pseudo-Riemannian manifold (implemented through Trac #24622). So if you upgrade to SageMath 8.2, you can declare an Euclidean space as a Riemannian manifold with a flat metric. It suffices to initialize the metric to the identity matrix in Cartesian coordinates:

```
sage: E = Manifold(3, 'E', structure='Riemannian')
sage: cartesian.<x,y,z> = E.chart()
sage: E.metric()[:] = identity_matrix(3)
```

Then, you can compute the curl of a vector field as follows:

```
sage: v = E.vector_field(name='v')
sage: v[:] = [-x*z, x+z, y^2]
sage: cv = v.curl()
sage: cv.display()
curl(v) = (2*y - 1) d/dx - x d/dy + d/dz
sage: cv[:]
[2*y - 1, -x, 1]
```

the divergence:

```
sage: dv = v.div()
sage: dv.display()
div(v): E --> R
(x, y, z) |--> -z
sage: dv.expr()
-z
```

the gradient:

```
sage: gdv = dv.gradient()
sage: gdv.display()
grad(div(v)) = -d/dz
sage: gdv[:]
[0, 0, -1]
```

the Laplacian:

```
sage: lv = v.laplacian()
sage: lv.display()
Delta(v) = 2 d/dz
sage: lv[:]
[0, 0, 2]
```

Instead of writing `v.div()`

, you may use the notation `div(v)`

:

```
sage: from sage.manifolds.operators import *
sage: div(v) == v.div()
True
```

Let us check a famous identity:

```
sage: curl(curl(v)) == grad(div(v)) - laplacian(v)
True
```