Motivation
Consider the following expression
$${\varepsilon}= \frac{1}{2} \left( \nabla \otimes u + \nabla \otimes u^\text{T} \right) \tag{1}$$
where $u:\mathbb{R^3} \to \mathbb{R^3}$ is a vector field, $\otimes$ is the tensor product, $\nabla$ is the the gradient operator and $\text{T}$ denotes the transpose. This expression rises in the theory of deformation of a continuum where it is assumed that the deformations are small. In fact, $u$ corresponds to a displacement field and $\varepsilon$ corresponds to the strain field which is a measure of deformation.
I want to compute the components of the strain field in cylindrical or spherical coordinates or any other curve-linear coordinates. So I go ahead by the usual computation rules in a curvelinear coordinates with basis $g_i$ and dual basis $g^i$. For example we can get
$$\begin{align} \nabla \otimes u &= (g^i \partial_i) \otimes (u_j g^j) \\ &= g^i \otimes \partial_i(u_j g^j) \\ &= g^i \otimes (\partial_i u_j g^j + u_j \partial_i g^j) \\ &= g^i \otimes (\partial_i u_j g^j -\Gamma_{ik}^{j} u_j g^k ) \\ &= g^i \otimes (\partial_i u_j g^j -\Gamma_{ij}^{k} u_k g^j ) \\ &= (\partial_i u_j -\Gamma_{ij}^{k} u_k ) g^i \otimes g^j \end{align} \tag{2}$$
Similarly, for the transpose we can get
$$ \nabla \otimes u ^ \text{T} = (\partial_j u_i -\Gamma_{ji}^{k} u_k ) g^i \otimes g^j \tag{3}$$
Finally, using the symmetry $\Gamma_{ij}^{k}=\Gamma_{ji}^{k}$ for Christoffel symbols we can get
$$\varepsilon=\frac{1}{2}(\partial_i u_j + \partial_j u_i -2 \Gamma_{ij}^{k} u_k ) g^i \otimes g^j \tag{4}$$
But this was really an easy case which we could do by hand. For example, consider the more complicated expression
$$\begin{align} t &= [\tau - (\nabla \cdot \mu)] \cdot n - \nabla \cdot (n \cdot \mu) \\ &+ [(\nabla \cdot n) (n \otimes n) + \nabla n^{\text{T}} \cdot (n \otimes n)]:\mu \\ &+ (n \otimes n \otimes n) \vdots \nabla \mu \end{align} \tag{5}$$
where $t$ and $n$ are vectors, $\tau$ is a second order tensor and $\mu$ is a third order tensor. The symbol $:$ is the double contraction and $\vdots$ is the triple contraction. It is more desirable to do this complex case with a CAS (Computer Algebra System) to avoid mistakes.
Questions
$1-$ How can I do these symbolic computations like in Eq. $(2)$ with a CAS like Mathematica or Maple?
$2-$ Is there already a code or software to support these kind of symbolic computations? If YES, can you kindly address it and do the symbolic computations in $(2)$ with the code or the software as an example in your answer.
While Mathematica is surely capable of handling abstract tensors/differential geometry computation, not so much capability is already built in but for special cases. Very recent versions have added tensor capabilities and the ability to define arbitrary coordinate systems, you may want to give it a try but I have no experience with it.
Another possibilities is to install some additional packages.
Atlas http://www.digi-area.com/Mathematica/atlas/guide/Atlas.php does the job. It is ok but, in my opinion, not great and it is not free. On the good side, it is not too complicated to use.
This one http://xact.es/xCoba/ looks great and is free, but seems pretty complicated to use and I have never got around using it.
Doubtless there are many other packages around. So bottom line is: for specific/sporadic applications you can probably get around with Mathematica built ins and a bit of programming. For more extensive use you are probably better of with a specific package (in Mathematica or any other CAS you like).