Is it possible to do vector calculus in Mathematica?

The answer depends a lot on what you mean by "doing" vector calculus. You want results to be displayed without using component notation, and that's in general a difficult thing to achieve.

A prerequisite about doing completely symbolic vector calculus is to define the simplification rules. But even in "non-vector" algebra it's often hard to get Simplify to reach an expression that you would like to obtain. This certainly won't get any easier when going to vector calculus. Trying to go the completely symbolic route may then not be worth your time, especially since I'm guessing you're doing everything in three dimensions, not "n dimensions."

So instead I usually follow a different approach in which Mathematica doesn't derive theorems but instead verifies them. That looks less impressive, but it's a lot easier to do. In your case, it means taking a vector calculus identity lhs == rhs and feeding it to Simplify in hopes of getting the result True. By contrast, starting only with Simplify[lhs], your chances of arriving exactly at rhs will in general by slim to nil.

To verify vector calculus identities, it's typically necessary to define your fields and coordinates in component form, but if you're lucky you won't have to display those components in the end result.

To show some examples, I wasn't able to make up my mind if I should use the VectorAnalysis package or the new version 9 functions. So I'll include my own versions of these functions here, just for the cartesian coordinate case.

Then I'll define something analogous to your vector r, and a general vector field notation. My choice for the coordinate symbols is to define them in a Array using a string "x" as the name for all components, so that there's no potential confusion with existing symbols.

Clear[grad, curl, div, cross, r];

grad[f_] := D[f, {Array["x", Length[f]]}]

div[f_] := Inner[D, f, Array["x", Length[f]]]

curl[f_?VectorQ] := 
 Inner[D, f.Normal[LeviCivitaTensor[Length[f]]], Array["x", Length[f]]]

cross[f_?VectorQ, g_?VectorQ] := 
 Dot[g.Normal[LeviCivitaTensor[Length[g]]], f]


r = Array["x", 3];

field[name_, dimension_: 3] := 
 Through[Array[name, dimension] @@ Array["x", dimension]]

Now check the vector identity that you mentioned in your comment:

Simplify[
 div[cross[field[A], field[B]]] == 
  field[B].curl[field[A]] - field[A].curl[field[B]]]

True

Next, check the continuity equation that motivated the question, by defining things in component form but displaying only the truth value of the identity:

Clear[jVec, ρ];

jVec[rVec_?VectorQ, t_] := 
 E^(-(rVec.rVec/(1 + t)^2))/(Pi^(3/2) (1 + t)^4) rVec

ρ[rVec_?VectorQ, t_] := 
 1/(Sqrt[Pi] Pi (1 + t)^3) Exp[-(rVec.rVec)/(1 + t)^2]

Simplify[
 div[jVec[r, t]] + D[ρ[r, t], t] == 0
 ]

True

When you do look at the coordinate form of the expressions, they will aways look somewhat uglier of course.

Our package VEST (Vector Einstein Summation Tools) is designed to do exactly this. It's described here and can be downloaded from github along with a tutorial.

I realize this question is now very old, but thought this answer might be helpful for others looking for similar capability.