Repository: sympy/sympy (python)

Issue:
Gradient of a vector

Per the discussion [here](https://groups.google.com/d/msgid/sympy/3f1832d4-b3b1-4241-b0be-6c02e7363359%40googlegroups.com?utm_medium=email&utm_source=footer), it appears as if the sympy `vector` module does not currently contain the ability to compute the gradient of a vector valued function. e.g. the Jacobian or Jacobian transpose. It would be nice if this capability existed since many PDEs involve vector variables, e.g. fluid/solid/continuum mechanics.
Improper Evaluation using sympy.vector.Del

When evaluating $\nabla \mathbf{v}$, where $\mathbf{v}$ is a vector (defined below), the `Del` operator appears to be incorrectly evaluated for cylindrical coordinates. (It is also incorrect in spherical coordinates).

```python
import sympy as sp
from sympy.vector import CoordSys3D, Del

R = CoordSys3D('R', transformation='cartesian', vector_names=("i", "j", "k")) # rectangular coordinates
C = CoordSys3D('C', transformation='cylindrical', vector_names=("e_r", "e_θ", "e_z"),
                                                  variable_names=('r', 'θ', 'z'))     # cylindrical coordinates

t = sp.symbols('t')                             # time

v_r = sp.Function('v_r')(t, C.r, C.θ, C.z)      # in general these components are functions of time and position
v_θ = sp.Function('v_θ')(t, C.r, C.θ, C.z)
v_z = sp.Function('v_z')(t, C.r, C.θ, C.z)

v = v_r*C.e_r + v_θ*C.e_θ + v_z*C.e_z

Del_v = Del()(v)

for coeff in [C.e_r, C.e_θ, C.e_z]:
    _ = Del()(v).coeff(coeff).expand()

    display(coeff)
    display(_)
    display(_.doit())
```

Output:

$$ \mathbf{\hat{e_r}_{C}} $$

$$ \left(\frac{\partial}{\partial \mathbf{{r}_{C}}} v_{r}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_r}_{C}} + \left(\frac{\partial}{\partial \mathbf{{r}_{C}}} v_{z}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_z}_{C}} + \left(\frac{\partial}{\partial \mathbf{{r}_{C}}} v_{\theta}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_\theta}_{C}} $$

$$ \mathbf{\hat{e_\theta}_{C}} $$

$$ \frac{\left(\frac{\partial}{\partial \mathbf{{\theta}_{C}}} v_{r}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_r}_{C}} + \left(\frac{\partial}{\partial \mathbf{{\theta}_{C}}} v_{z}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_z}_{C}} + \left(\frac{\partial}{\partial \mathbf{{\theta}_{C}}} v_{\theta}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_\theta}_{C}}}{\mathbf{{r}_{C}}} $$

$$ \mathbf{\hat{e_z}_{C}} $$

$$ \left(\frac{\partial}{\partial \mathbf{{z}_{C}}} v_{r}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_r}_{C}} + \left(\frac{\partial}{\partial \mathbf{{z}_{C}}} v_{z}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_z}_{C}} + \left(\frac{\partial}{\partial \mathbf{{z}_{C}}} v_{\theta}{\left(t,\mathbf{{r}_{C}},\mathbf{{\theta}_{C}},\mathbf{{z}_{C}} \right)}\right)\mathbf{\hat{e_\theta}_{C}} $$

The output above is correct in the $r$ and $z$ directions, but in the $\theta$ direction, the derivatives of the unit vectors 
$ \mathbf{e}_r $ and $ \mathbf{e}_\theta $ are not properly evaluated. 

$$ \mathbf{e}_r = \cos{(\theta)} \mathbf{e}_x + \sin{(\theta)} \mathbf{e}_y $$

$$ \mathbf{e}_\theta = -\sin{(\theta)} \mathbf{e}_x + \cos{(\theta)} \mathbf{e}_y $$

$$ \frac{\partial}{\partial \theta} \left(\mathbf{e}_r \right) = -\sin{(\theta)} \mathbf{e}_x + \cos{(\theta)} \mathbf{e}_y = \mathbf{e}_\theta $$

$$ \frac{\partial}{\partial \theta} \left(\mathbf{e}_\theta \right) = -\cos{(\theta)} \mathbf{e}_x - \sin{(\theta)} \mathbf{e}_y = -\mathbf{e}_r $$

Where as `Del()(v)` seems to evaluate the derivatives as: 

$$ \frac{\partial}{\partial \theta} \left(\mathbf{e}_r \right) = \frac{\partial}{\partial \theta} \left(\mathbf{e}_\theta \right)  = 0 $$
field-aligned gradient of vectorfield does not give correct results in cylindrical coordinate systems

To calculate terms in fluid dynamical differential equations in the style of  $\vec{u}\cdot\nabla\vec{v}$, none of the available implementations in sympy give the correct result for a cylindrical coordinate system. I have tried the two implementations as suggested in issue #12969 and issue #19212.
In the example below, I use the vector field $\vec{v} = (0,\Omega,0)$ in cylindrical coordinates. The expected result for $\vec{v}\cdot\nabla\vec{v}$ is $(-\Omega^2/r,0,0)$ (including the pesky $-A_\phi B_\phi/r$ terms as can be found on p. 7 of [the NRL plasma formulary](https://library.psfc.mit.edu/catalog/online_pubs/NRL_FORMULARY_13.pdf)). 

```python
import sympy as sp 
from sympy.vector import CoordSys3D, gradient, curl, divergence, directional_derivative
c = CoordSys3D('c', transformation='cylindrical', variable_names=("r", "th", "z"))
Omega = sp.symbols('Omega')
test = Omega*c.j
test2=directional_derivative(test, test) # this does not yield the correct expression: it should be (-Omega^2/r, 0, 0)
test3=test&gradient(test) # this does not yield the correct expression: it even remains a scalar, even though it should be a vector 
```
``test2`` is just the $\vec{0}$ vector, and ``test3`` is a 0 scalar. Neither are correct. 

A workaround is possible using the ``galgebra`` package. Here is a working example
```python
from galgebra.ga import Ga
X = (r,th,z) = sp.symbols('r theta z')
c3d = Ga('e_r e_theta e_z',g=[1,r**2,1],coords=X,norm=True)
(er,eth,ez) = c3d.mv()
grad = c3d.grad
test = Omega*eth
test4=(test|grad)*test
```
which does return the expected result. 

PS. I think the same remark holds for spherical coordinates. 


Patch:
diff --git a/sympy/vector/__init__.py b/sympy/vector/__init__.py
index f6757bbeb350..d3e1c120edc7 100644
--- a/sympy/vector/__init__.py
+++ b/sympy/vector/__init__.py
@@ -6,6 +6,7 @@
 from sympy.vector.scalar import BaseScalar
 from sympy.vector.deloperator import Del
 from sympy.vector.functions import (express, matrix_to_vector,
+                                    matrix_to_dyadic,
                                     laplacian, is_conservative,
                                     is_solenoidal, scalar_potential,
                                     directional_derivative,
@@ -33,7 +34,8 @@
 
     'CoordSys3D',
 
-    'express', 'matrix_to_vector', 'laplacian', 'is_conservative',
+    'express', 'matrix_to_vector', 'matrix_to_dyadic',
+    'laplacian', 'is_conservative',
     'is_solenoidal', 'scalar_potential', 'directional_derivative',
     'scalar_potential_difference',
 
diff --git a/sympy/vector/functions.py b/sympy/vector/functions.py
index 252c976ba797..cc6105e174c0 100644
--- a/sympy/vector/functions.py
+++ b/sympy/vector/functions.py
@@ -1,13 +1,14 @@
+from sympy.matrices.exceptions import ShapeError
 from sympy.vector.coordsysrect import CoordSys3D
 from sympy.vector.deloperator import Del
 from sympy.vector.scalar import BaseScalar
-from sympy.vector.vector import Vector, BaseVector
-from sympy.vector.operators import gradient, curl, divergence
+from sympy.vector.vector import Vector, BaseVector, Dot
+from sympy.vector.operators import gradient, curl, divergence, Gradient
 from sympy.core.function import diff
 from sympy.core.singleton import S
 from sympy.integrals.integrals import integrate
 from sympy.core import sympify
-from sympy.vector.dyadic import Dyadic
+from sympy.vector.dyadic import Dyadic, DyadicAdd
 
 
 def express(expr, system, system2=None, variables=False):
@@ -125,9 +126,34 @@ def express(expr, system, system2=None, variables=False):
 
 
 def directional_derivative(field, direction_vector):
-    """
+    r"""
     Returns the directional derivative of a scalar or vector field computed
-    along a given vector in coordinate system which parameters are expressed.
+    along a given vector in a coordinate system.
+
+    Note that there are two definitions of directional derivative for a
+    differentiable scalar function $f$. Let $\mathbf{\vec{v}}$ be an arbitrary
+    direction vector, and $\mathbf{\hat{v}}$ be the normalized direction
+    vector. Some authors define the directional derivative as:
+
+    $$\nabla_{\mathbf{\vec{v}}} f(\mathbf{x}) = \mathbf{\hat{v}} \cdot \nabla f(\mathbf{x})$$
+
+    Other authors define it as:
+
+    $$\nabla_{\vec{v}} f(\mathbf{x}) = \vec{v} \cdot \nabla f(\mathbf{x})$$
+
+    The two definitions differs by the magnitude of the direction vector.
+    This function does not normalize direction_vector. If the user wants
+    to compute the directional derivative according to the first definition,
+    it is the user's resposibility to normalize the direction vector before
+    giving it to this function.
+
+    The directional derivative of a vector field $\mathbf{\vec{f}}(\mathbf{x})$
+    along the vector $\mathbf{\vec{v}}$ is:
+
+    $$\nabla_{\mathbf{\vec{v}}} f(\mathbf{x}) = \mathbf{\vec{v}} \cdot \nabla \mathbf{\vec{f}}(\mathbf{x})$$
+
+    This term frequently appears in continuum mechanics's partial differential
+    equations, in particular the advection equation.
 
     Parameters
     ==========
@@ -138,40 +164,55 @@ def directional_derivative(field, direction_vector):
     direction_vector : Vector
         The vector to calculated directional derivative along them.
 
-
     Examples
     ========
 
+    >>> from sympy import Function
     >>> from sympy.vector import CoordSys3D, directional_derivative
     >>> R = CoordSys3D('R')
+
+    Directional derivative of a scalar field with an arbitrary direction
+    vector:
+
     >>> f1 = R.x*R.y*R.z
     >>> v1 = 3*R.i + 4*R.j + R.k
     >>> directional_derivative(f1, v1)
     R.x*R.y + 4*R.x*R.z + 3*R.y*R.z
-    >>> f2 = 5*R.x**2*R.z
-    >>> directional_derivative(f2, v1)
-    5*R.x**2 + 30*R.x*R.z
+
+    Directional derivative of a scalar field with a unit vector:
+
+    >>> directional_derivative(f1, v1.normalize())
+    sqrt(26)*R.x*R.y/26 + 2*sqrt(26)*R.x*R.z/13 + 3*sqrt(26)*R.y*R.z/26
+
+    Directional derivate of a vector field in a cylindrical coordinate system:
+
+    >>> C = CoordSys3D("C", transformation="cylindrical")
+    >>> r, theta, z = C.base_scalars()
+    >>> e_r, e_theta, e_z = C.base_vectors()
+    >>> u, v, w = [Function(s)(r, theta, z) for s in ["u", "v", "w"]]
+    >>> vec = u * e_r + v * e_theta + w * e_z
+    >>> direction = e_theta
+    >>> directional_derivative(vec, direction)  # doctest: +NORMALIZE_WHITESPACE
+    (-v(C.r, C.theta, C.z)/C.r + Derivative(u(C.r, C.theta, C.z), C.theta)/C.r)*C.i
+    + (u(C.r, C.theta, C.z)/C.r + Derivative(v(C.r, C.theta, C.z), C.theta)/C.r)*C.j
+    + (Derivative(w(C.r, C.theta, C.z), C.theta)/C.r)*C.k
 
     """
     from sympy.vector.operators import _get_coord_systems
     coord_sys = _get_coord_systems(field)
-    if len(coord_sys) > 0:
-        # TODO: This gets a random coordinate system in case of multiple ones:
-        coord_sys = next(iter(coord_sys))
-        field = express(field, coord_sys, variables=True)
-        i, j, k = coord_sys.base_vectors()
-        x, y, z = coord_sys.base_scalars()
-        out = Vector.dot(direction_vector, i) * diff(field, x)
-        out += Vector.dot(direction_vector, j) * diff(field, y)
-        out += Vector.dot(direction_vector, k) * diff(field, z)
-        if out == 0 and isinstance(field, Vector):
-            out = Vector.zero
-        return out
-    elif isinstance(field, Vector):
-        return Vector.zero
-    else:
+    if len(coord_sys) == 0:
+        if isinstance(field, Vector):
+            return Vector.zero
         return S.Zero
 
+    coord_sys = next(iter(coord_sys))
+    field = express(field, coord_sys, variables=True)
+    grad = gradient(field)
+
+    if isinstance(grad, Gradient):
+        return Dot(direction_vector, grad)
+    return direction_vector & grad
+
 
 def laplacian(expr):
     """
@@ -429,6 +470,46 @@ def matrix_to_vector(matrix, system):
     return outvec
 
 
+def matrix_to_dyadic(matrix, system):
+    """
+    Converts a 3 x 3 matrix to a Dyadic instance.
+
+    Parameters
+    ==========
+
+    matrix : Matrix, Dimensions: (3, 3)
+        The matrix to be converted to a dyadic
+
+    system : CoordSys3D
+        The coordinate system the dyadic is to be defined in
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> from sympy.vector import CoordSys3D, matrix_to_dyadic
+    >>> C = CoordSys3D('C')
+    >>> m = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+    >>> d = matrix_to_dyadic(m, C)
+    >>> d     # doctest: +NORMALIZE_WHITESPACE
+    (C.i|C.i) + 2*(C.i|C.j) + 3*(C.i|C.k) + 4*(C.j|C.i) + 5*(C.j|C.j)
+    + 6*(C.j|C.k) + 7*(C.k|C.i) + 8*(C.k|C.j) + 9*(C.k|C.k)
+    >>> d.to_matrix(C) == m
+    True
+
+    """
+    n, m = matrix.shape
+    if (n != 3) or (m != 3):
+        raise ShapeError("A 3 x 3 matrix is required.")
+
+    args = []
+    vects = system.base_vectors()
+    for i in range(3):
+        for j in range(3):
+            args.append(matrix[i, j] * (vects[i] | vects[j]))
+    return DyadicAdd(*args)
+
+
 def _path(from_object, to_object):
     """
     Calculates the 'path' of objects starting from 'from_object'
diff --git a/sympy/vector/operators.py b/sympy/vector/operators.py
index 3ca42d20302f..e2daa431c08b 100644
--- a/sympy/vector/operators.py
+++ b/sympy/vector/operators.py
@@ -2,6 +2,7 @@
 from sympy.core.expr import Expr
 from sympy.core import sympify, S, preorder_traversal
 from sympy.vector.coordsysrect import CoordSys3D
+from sympy.vector.dyadic import Dyadic, DyadicAdd
 from sympy.vector.vector import Vector, VectorMul, VectorAdd, Cross, Dot
 from sympy.core.function import Derivative
 from sympy.core.add import Add
@@ -242,16 +243,32 @@ def divergence(vect, doit=True):
             raise ValueError("Invalid argument for divergence")
 
 
-def gradient(scalar_field, doit=True):
+def _christoffel_symbol_2nd_kind(h, q, i, j, k):
     """
-    Returns the vector gradient of a scalar field computed wrt the
-    base scalars of the given coordinate system.
+    Compute the Christoffel symbol of the second kind, Gamma^{k}_{ij},
+    for an orthogonal system.
+    """
+    if (i == j) and (i == k):
+        return Derivative(h[i], q[i]) / h[i]
+    if (i == k) and (i != j):
+        return Derivative(h[i], q[j]) / h[i]
+    if (j == k) and (i != j):
+        return Derivative(h[j], q[i]) / h[j]
+    if (i == j) and (i != k):
+        return -Derivative(h[i], q[k]) * h[i] / h[k]**2
+    return S.Zero
+
+
+def gradient(field, doit=True):
+    """
+    Returns the gradient of the field computed wrt the base scalars
+    of the given coordinate system.
 
     Parameters
     ==========
 
-    scalar_field : SymPy Expr
-        The scalar field to compute the gradient of
+    field : Expr, Vector
+        The scalar-valued or vector-valued function representing the field.
 
     doit : bool
         If True, the result is returned after calling .doit() on
@@ -261,7 +278,11 @@ def gradient(scalar_field, doit=True):
     Examples
     ========
 
+    >>> from sympy import Function
     >>> from sympy.vector import CoordSys3D, gradient
+
+    The gradient of a scalar field is a vector field:
+
     >>> R = CoordSys3D('R')
     >>> s1 = R.x*R.y*R.z
     >>> gradient(s1)
@@ -270,30 +291,85 @@ def gradient(scalar_field, doit=True):
     >>> gradient(s2)
     10*R.x*R.z*R.i + 5*R.x**2*R.k
 
-    """
-    coord_sys = _get_coord_systems(scalar_field)
+    The gradient of a vector field is a dyadic (second order tensor):
+
+    >>> C = CoordSys3D("C", transformation="cylindrical")
+    >>> e_r, e_theta, e_z = C.base_vectors()
+    >>> r, theta, z = C.base_scalars()
+    >>> v_r = Function('v_r')(r, theta, z)
+    >>> v_theta = Function('v_theta')(r, theta, z)
+    >>> v_z = Function('v_z')(r, theta, z)
+    >>> v = v_r * e_r + v_theta * e_theta + v_z * e_z
+    >>> gradient(v)     # doctest: +NORMALIZE_WHITESPACE
+    (Derivative(v_r(C.r, C.theta, C.z), C.r))*(C.i|C.i)
+    + (Derivative(v_theta(C.r, C.theta, C.z), C.r))*(C.i|C.j)
+    + (Derivative(v_z(C.r, C.theta, C.z), C.r))*(C.i|C.k)
+    + (-v_theta(C.r, C.theta, C.z)/C.r + Derivative(v_r(C.r, C.theta, C.z), C.theta)/C.r)*(C.j|C.i)
+    + (v_r(C.r, C.theta, C.z)/C.r + Derivative(v_theta(C.r, C.theta, C.z), C.theta)/C.r)*(C.j|C.j)
+    + (Derivative(v_z(C.r, C.theta, C.z), C.theta)/C.r)*(C.j|C.k)
+    + (Derivative(v_r(C.r, C.theta, C.z), C.z))*(C.k|C.i)
+    + (Derivative(v_theta(C.r, C.theta, C.z), C.z))*(C.k|C.j)
+    + (Derivative(v_z(C.r, C.theta, C.z), C.z))*(C.k|C.k)
 
+    """
+    coord_sys = _get_coord_systems(field)
     if len(coord_sys) == 0:
         return Vector.zero
-    elif len(coord_sys) == 1:
+
+    if not isinstance(field, (Vector, Dyadic)):
+        # scalar field
+        if len(coord_sys) > 1:
+            if isinstance(field, Add):
+                return VectorAdd.fromiter(gradient(i) for i in field.args)
+            if isinstance(field, Mul):
+                s = _split_mul_args_wrt_coordsys(field)
+                return VectorAdd.fromiter(field / i * gradient(i) for i in s)
+            return Gradient(field)
+
         coord_sys = next(iter(coord_sys))
         h1, h2, h3 = coord_sys.lame_coefficients()
         i, j, k = coord_sys.base_vectors()
         x, y, z = coord_sys.base_scalars()
-        vx = Derivative(scalar_field, x) / h1
-        vy = Derivative(scalar_field, y) / h2
-        vz = Derivative(scalar_field, z) / h3
+        vx = Derivative(field, x) / h1
+        vy = Derivative(field, y) / h2
+        vz = Derivative(field, z) / h3
 
         if doit:
             return (vx * i + vy * j + vz * k).doit()
         return vx * i + vy * j + vz * k
-    else:
-        if isinstance(scalar_field, (Add, VectorAdd)):
-            return VectorAdd.fromiter(gradient(i) for i in scalar_field.args)
-        if isinstance(scalar_field, (Mul, VectorMul)):
-            s = _split_mul_args_wrt_coordsys(scalar_field)
-            return VectorAdd.fromiter(scalar_field / i * gradient(i) for i in s)
-        return Gradient(scalar_field)
+
+    elif isinstance(field, Vector):
+        # vector field
+        if len(coord_sys) > 1:
+            if isinstance(field, VectorAdd):
+                return DyadicAdd.fromiter(gradient(i) for i in field.args)
+            if isinstance(field, VectorMul):
+                s = _split_mul_args_wrt_coordsys(field)
+                return DyadicAdd.fromiter(field / i * gradient(i) for i in s)
+            return Gradient(field)
+
+        coord_sys = next(iter(coord_sys))
+        e = coord_sys.base_vectors()
+        q = coord_sys.base_scalars()
+        h = coord_sys.lame_coefficients()
+        A = [field & u for u in e]
+        s = Dyadic.zero
+
+        for i in range(3):
+            for j in range(3):
+                component = Derivative(A[j], q[i]) / h[i] - \
+                    A[j] / (h[i] * h[j]) * Derivative(h[j], q[i])
+
+                for k in range(3):
+                    Gamma = _christoffel_symbol_2nd_kind(h, q, i, k, j)
+                    component += h[j] / (h[i] * h[k]) * A[k] * Gamma
+
+                s += component * (e[i] | e[j])
+        if doit:
+            return s.doit()
+        return s
+
+    return Gradient(field)
 
 
 class Laplacian(Expr):


Test patch:
diff --git a/sympy/testing/quality_unicode.py b/sympy/testing/quality_unicode.py
index f8d7d888a85c..ed2dfcbdf2e5 100644
--- a/sympy/testing/quality_unicode.py
+++ b/sympy/testing/quality_unicode.py
@@ -27,6 +27,7 @@
     r'*/sympy/physics/vector/tests/test_printing.py',
     r'*/physics/quantum/tests/test_printing.py',
     r'*/sympy/vector/tests/test_printing.py',
+    r'*/sympy/vector/tests/test_field_functions.py',
     r'*/sympy/parsing/tests/test_sympy_parser.py',
     r'*/sympy/printing/pretty/stringpict.py',
     r'*/sympy/printing/pretty/tests/test_pretty.py',
diff --git a/sympy/vector/tests/test_field_functions.py b/sympy/vector/tests/test_field_functions.py
index 035c2ce0234b..4bfabef440bb 100644
--- a/sympy/vector/tests/test_field_functions.py
+++ b/sympy/vector/tests/test_field_functions.py
@@ -1,4 +1,6 @@
-from sympy.core.function import Derivative
+from sympy.core.function import Derivative, Function, diff
+from sympy.matrices import Matrix
+from sympy.vector.dyadic import Dyadic
 from sympy.vector.vector import Vector
 from sympy.vector.coordsysrect import CoordSys3D
 from sympy.simplify import simplify
@@ -8,10 +10,13 @@
 from sympy.vector.vector import Dot
 from sympy.vector.operators import curl, divergence, gradient, Gradient, Divergence, Cross
 from sympy.vector.deloperator import Del
-from sympy.vector.functions import (is_conservative, is_solenoidal,
-                                    scalar_potential, directional_derivative,
-                                    laplacian, scalar_potential_difference)
+from sympy.vector.functions import (
+        is_conservative, is_solenoidal,
+        scalar_potential, directional_derivative,
+        laplacian, scalar_potential_difference,
+        matrix_to_dyadic)
 from sympy.testing.pytest import raises
+import functools
 
 C = CoordSys3D('C')
 i, j, k = C.base_vectors()
@@ -20,6 +25,36 @@
 a, b, c, q = symbols('a b c q')
 
 
+@functools.cache
+def _setup_cartesian_system():
+    C = CoordSys3D("C")
+    x, y, z = C.base_scalars()
+    i, j, k = C.base_vectors()
+    u, v, w = [Function(s)(x, y, z) for s in ["u", "v", "w"]]
+    vec = u * i + v * j + w * k
+    return C, x, y, z, i, j, k, u, v, w, vec
+
+
+@functools.cache
+def _setup_cylindrical_system():
+    C = CoordSys3D("C", transformation="cylindrical")
+    r, θ, z = C.base_scalars()
+    e_r, e_θ, e_z = C.base_vectors()
+    u, v, w = [Function(s)(r, θ, z) for s in ["u", "v", "w"]]
+    vec = u * e_r + v * e_θ + w * e_z
+    return C, r, θ, z, e_r, e_θ, e_z, u, v, w, vec
+
+
+@functools.cache
+def _setup_spherical_system():
+    S = CoordSys3D("S", transformation="spherical")
+    r, θ, 𐌘 = S.base_scalars()
+    e_r, e_θ, e_𐌘 = S.base_vectors()
+    u, v, w = [Function(s)(r, θ, 𐌘) for s in ["u", "v", "w"]]
+    vec = u * e_r + v * e_θ + w * e_𐌘
+    return S, r, θ, 𐌘, e_r, e_θ, e_𐌘, u, v, w, vec
+
+
 def test_del_operator():
     # Tests for curl
 
@@ -211,7 +246,49 @@ def test_directional_derivative():
     r, theta, phi = D.base_scalars()
     e_r, e_theta, e_phi = D.base_vectors()
     assert directional_derivative(r**2*e_r, e_r) == 2*r*e_r
-    assert directional_derivative(5*r**2*phi, 3*e_r + 4*e_theta + e_phi) == 5*r**2 + 30*r*phi
+    assert directional_derivative(5*r**2*phi, 3*e_r + 4*e_theta + e_phi) == 5*r/sin(theta) + 30*r*phi
+
+
+def test_directional_derivative_vector_field_cartesian_systems():
+    C, x, y, z, i, j, k, u, v, w, vec = _setup_cartesian_system()
+    ux, uy, uz = [u.diff(s) for s in [x, y, z]]
+    vx, vy, vz = [v.diff(s) for s in [x, y, z]]
+    wx, wy, wz = [w.diff(s) for s in [x, y, z]]
+    direction = i + j + k
+
+    res = directional_derivative(vec, direction)
+    assert res == (
+        (ux + uy + uz) * i +
+        (vx + vy + vz) * j +
+        (wx + wy + wz) * k)
+
+
+def test_directional_derivative_vector_field_cylindrical_systems():
+    C, r, θ, z, e_r, e_θ, e_z, u, v, w, vec = _setup_cylindrical_system()
+    ur, ut, uz = [u.diff(s) for s in [r, θ, z]]
+    vr, vt, vz = [v.diff(s) for s in [r, θ, z]]
+    wr, wt, wz = [w.diff(s) for s in [r, θ, z]]
+    direction = e_r + e_θ + e_z
+
+    res = directional_derivative(vec, direction)
+    assert res == (
+        (ur + ut / r + uz - v / r) * e_r +
+        (vr + vt / r + vz + u / r) * e_θ +
+        (wr + wt / r + wz) * e_z)
+
+
+def test_directional_derivative_vector_field_spherical_systems():
+    S, r, θ, 𐌘, e_r, e_θ, e_𐌘, u, v, w, vec = _setup_spherical_system()
+    ur, ut, up = [u.diff(s) for s in [r, θ, 𐌘]]
+    vr, vt, vp = [v.diff(s) for s in [r, θ, 𐌘]]
+    wr, wt, wp = [w.diff(s) for s in [r, θ, 𐌘]]
+    direction = e_r + e_θ + e_𐌘
+
+    res = directional_derivative(vec, direction)
+    assert res == (
+        (ur + ut / r + up / (r * sin(θ)) - v / r - w / r) * e_r +
+        (vr + vt / r + vp / (r * sin(θ)) + u / r - w * cos(θ) / (r * sin(θ))) * e_θ +
+        (wr + wt / r + wp / (r * sin(θ)) + u / r + v * cos(θ) / (r * sin(θ))) * e_𐌘)
 
 
 def test_scalar_potential():
@@ -319,3 +396,187 @@ def test_mixed_coordinates():
                 a.x**2*b.x*Dot(a.i, c.i) +\
                 b.x**2*c.x*Dot(b.i, a.i) +\
                 a.x*b.x**2*Dot(b.i, c.i)
+
+
+def test_gradient_of_vector_cartesian_coordinates():
+    C, x, y, z, i, j, k, u, v, w, vec = _setup_cartesian_system()
+
+    res = gradient(vec)
+    assert isinstance(res, Dyadic)
+    assert res.to_matrix(C) == Matrix([
+        [diff(u, x), diff(v, x), diff(w, x)],
+        [diff(u, y), diff(v, y), diff(w, y)],
+        [diff(u, z), diff(v, z), diff(w, z)]
+    ])
+    assert delop(vec).doit() == res
+
+
+def test_gradient_of_vector_cylindrical_coordinates():
+    C, r, θ, z, e_r, e_θ, e_z, u, v, w, vec = _setup_cylindrical_system()
+
+    res = gradient(vec)
+    assert isinstance(res, Dyadic)
+    assert res.to_matrix(C) == Matrix([
+        [diff(u, r), diff(v, r), diff(w, r)],
+        [diff(u, θ) / r - v / r, diff(v, θ) / r + u / r, diff(w, θ) / r],
+        [diff(u, z), diff(v, z), diff(w, z)]
+    ])
+    assert delop(vec).doit() == res
+
+
+def test_gradient_of_vector_spherical_coordinates():
+    S, r, θ, 𐌘, e_r, e_θ, e_𐌘, u, v, w, vec = _setup_spherical_system()
+
+    res = gradient(vec)
+    assert isinstance(res, Dyadic)
+    assert res.to_matrix(S) == Matrix([
+        [diff(u, r), diff(v, r), diff(w, r)],
+        [diff(u, θ) / r - v / r, diff(v, θ) / r + u / r, diff(w, θ) / r],
+        [
+                diff(u, 𐌘) / (r * sin(θ)) - w / r,
+                diff(v, 𐌘) / (r * sin(θ)) - w * cos(θ)  / (r * sin(θ)),
+                u / r + v * cos(θ) / (r * sin(θ)) + diff(w, 𐌘) / (r * sin(θ))
+        ]
+    ])
+    assert delop(vec).doit() == res
+
+
+def test_gradient_of_vector_cartesian_cylindrical():
+    # verify that gradients works as expected for a vector whose
+    # components are defined in multiple coordinate systems
+
+    Cart = CoordSys3D("Cart")
+    x, y, z = Cart.base_scalars()
+    v_x = Function('v_x')(x, y, z)
+    v_y = Function('v_y')(x, y, z)
+    v_z = Function('v_z')(x, y, z)
+    v1 = v_x * Cart.i + v_y * Cart.j + v_z * Cart.k
+    m1 = Matrix([
+        [diff(v_x, x), diff(v_y, x), diff(v_z, x)],
+        [diff(v_x, y), diff(v_y, y), diff(v_z, y)],
+        [diff(v_x, z), diff(v_y, z), diff(v_z, z)]
+    ])
+    d1 = matrix_to_dyadic(m1, Cart)
+
+    C = Cart.create_new("C", transformation="cylindrical")
+    r, θ, z = C.base_scalars()
+    v_r = Function('v_r')(r, θ, z)
+    v_θ = Function('v_θ')(r, θ, z)
+    v_z = Function('v_z')(r, θ, z)
+    v2 = v_r * C.i + v_θ * C.j + v_z * C.k
+    m2 = Matrix([
+        [diff(v_r, r), diff(v_θ, r), diff(v_z, r)],
+        [diff(v_r, θ) / r - v_θ / r, diff(v_θ, θ) / r + v_r / r, diff(v_z, θ) / r],
+        [diff(v_r, z), diff(v_θ, z), diff(v_z, z)]
+    ])
+    d2 = matrix_to_dyadic(m2, C)
+
+    v = v1 + v2
+    res = gradient(v)
+    assert res == d1 + d2
+    assert delop(v).doit() == res
+
+
+def test_issue_27427():
+    C = CoordSys3D('C', transformation='cylindrical')
+    r, theta, z = C.base_scalars()
+    Omega = symbols('Omega')
+    v = Omega * C.j
+    expected = -Omega**2 / r * C.i
+    assert v & gradient(v) == expected
+    assert directional_derivative(v, v) == expected
+
+
+def test_laplacian_vector_field_cartesian_system():
+    C, x, y, z, i, j, k, u, v, w, vec = _setup_cartesian_system()
+
+    assert laplacian(vec) == (
+        (u.diff(x, 2) + u.diff(y, 2) + u.diff(z, 2)) * i +
+        (v.diff(x, 2) + v.diff(y, 2) + v.diff(z, 2)) * j +
+        (w.diff(x, 2) + w.diff(y, 2) + w.diff(z, 2)) * k)
+
+
+def test_laplacian_vector_field_cylindrical_system():
+    C, r, θ, z, e_r, e_θ, e_z, u, v, w, vec = _setup_cylindrical_system()
+
+    assert laplacian(vec).expand() == (
+        (u.diff(r, 2) + u.diff(θ, 2) / r**2 + u.diff(z, 2) + u.diff(r) / r - v.diff(θ) * 2 / r**2 - u / r**2) * e_r +
+        (v.diff(r, 2) + v.diff(θ, 2) / r**2 + v.diff(z, 2) + v.diff(r) / r + u.diff(θ) * 2 / r**2 - v / r**2) * e_θ +
+        (w.diff(r, 2) + w.diff(θ, 2) / r**2 + w.diff(z, 2) + w.diff(r) / r) * e_z)
+
+
+def test_laplacian_vector_field_spherical_system():
+    S, r, θ, 𐌘, e_r, e_θ, e_𐌘, u, v, w, vec = _setup_spherical_system()
+
+    expected_c1 = (
+        ((r**2 * u).diff(r) / r**2).diff(r)
+        + (sin(θ) * u.diff(θ)).diff(θ) / (r**2 * sin(θ))
+        + u.diff(𐌘, 2) / (r**2 * sin(θ)**2)
+        - 2 * (v * sin(θ)).diff(θ) / (r**2 * sin(θ))
+        - 2 * w.diff(𐌘) / (r**2 * sin(θ))).expand()
+    expected_c2 = (
+        (r**2 * v.diff(r)).diff(r) / r**2
+        + ((v * sin(θ)).diff(θ) / sin(θ)).diff(θ) / r**2
+        + v.diff(𐌘, 2) / (r**2 * sin(θ)**2)
+        + 2 * u.diff(θ) / r**2
+        - 2 * cos(θ) / (r**2 * sin(θ)**2) * w.diff(𐌘)).expand()
+    expected_c3 = (
+        (w.diff(r) * r**2).diff(r) / r**2
+        + ((w * sin(θ)).diff(θ) / sin(θ)).diff(θ) / r**2
+        + w.diff(𐌘, 2) / (r**2 * sin(θ)**2)
+        + 2 * u.diff(𐌘) / (r**2 * sin(θ))
+        + 2 * cos(θ) / (r**2 * sin(θ)**2) * v.diff(𐌘)).expand()
+    assert laplacian(vec).expand() == (
+        expected_c1 * e_r +
+        expected_c2 * e_θ +
+        expected_c3 * e_𐌘)
+
+
+def test_divergence_cartesian_system():
+    C, x, y, z, i, j, k, u, v, w, vec = _setup_cartesian_system()
+
+    assert divergence(vec) == u.diff(x) + v.diff(y) + w.diff(z)
+
+
+def test_divergence_cylindrical_system():
+    C, r, θ, z, e_r, e_θ, e_z, u, v, w, vec = _setup_cylindrical_system()
+
+    expected = ((r * u).diff(r) + v.diff(θ) + r * w.diff(z)) / r
+    assert divergence(vec).expand() == expected.expand()
+
+
+def test_divergence_spherical_system():
+    S, r, θ, 𐌘, e_r, e_θ, e_𐌘, u, v, w, vec = _setup_spherical_system()
+
+    expected = (
+        sin(θ) * (r**2 * u).diff(r)
+        + r * (v * sin(θ)).diff(θ)
+        + r * w.diff(𐌘)) / (r**2 * sin(θ))
+    assert divergence(vec).expand() == expected.expand()
+
+
+def test_curl_cartesian_system():
+    C, x, y, z, i, j, k, u, v, w, vec = _setup_cartesian_system()
+
+    assert curl(vec) == (
+        (w.diff(y) - v.diff(z)) * i +
+        (u.diff(z) - w.diff(x)) * j +
+        (v.diff(x) - u.diff(y)) * k)
+
+
+def test_curl_cylindrical_system():
+    C, r, θ, z, e_r, e_θ, e_z, u, v, w, vec = _setup_cylindrical_system()
+
+    assert curl(vec).expand() == (
+        (w.diff(θ) / r - v.diff(z)) * e_r +
+        (u.diff(z) - w.diff(r)) * e_θ +
+        (v / r + v.diff(r) - u.diff(θ) / r) * e_z)
+
+
+def test_curl_spherical_system():
+    S, r, θ, 𐌘, e_r, e_θ, e_𐌘, u, v, w, vec = _setup_spherical_system()
+
+    assert curl(vec).expand() == (
+        (w * cos(θ) / (r * sin(θ)) + w.diff(θ) / r - v.diff(𐌘) / (r * sin(θ))) * e_r +
+        (u.diff(𐌘) / (r * sin(θ)) - w / r - w.diff(r)) * e_θ +
+        (v.diff(r) + v / r - u.diff(θ) / r) * e_𐌘)
diff --git a/sympy/vector/tests/test_functions.py b/sympy/vector/tests/test_functions.py
index 87fb871b52fc..50d893c0c517 100644
--- a/sympy/vector/tests/test_functions.py
+++ b/sympy/vector/tests/test_functions.py
@@ -1,11 +1,13 @@
 from sympy.vector.vector import Vector
 from sympy.vector.coordsysrect import CoordSys3D
-from sympy.vector.functions import express, matrix_to_vector, orthogonalize
+from sympy.vector.functions import (
+    express, matrix_to_vector, matrix_to_dyadic, orthogonalize)
 from sympy.core.numbers import Rational
 from sympy.core.singleton import S
 from sympy.core.symbol import symbols
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import (cos, sin)
+from sympy.matrices.exceptions import ShapeError
 from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
 from sympy.testing.pytest import raises
 
@@ -165,6 +167,21 @@ def test_matrix_to_vector():
     assert matrix_to_vector(m, N) == q1*N.i + q2*N.j + q3*N.k
 
 
+def test_matrix_to_dyadic():
+    raises(ShapeError, lambda: matrix_to_dyadic(Matrix([0]), N))
+    raises(
+        ShapeError,
+        lambda: matrix_to_dyadic(Matrix([[0, 1, 2], [3, 4, 5]]), N))
+
+    m = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+    expected = (
+        1 * (N.i | N.i) + 2 * (N.i | N.j) + 3 * (N.i | N.k) +
+        4 * (N.j | N.i) + 5 * (N.j | N.j) + 6 * (N.j | N.k) +
+        7 * (N.k | N.i) + 8 * (N.k | N.j) + 9 * (N.k | N.k)
+    )
+    assert matrix_to_dyadic(m, N) == expected
+
+
 def test_orthogonalize():
     C = CoordSys3D('C')
     a, b = symbols('a b', integer=True)
diff --git a/sympy/vector/tests/test_operators.py b/sympy/vector/tests/test_operators.py
index 5734edadd005..92749730abeb 100644
--- a/sympy/vector/tests/test_operators.py
+++ b/sympy/vector/tests/test_operators.py
@@ -1,5 +1,8 @@
 from sympy.vector import CoordSys3D, Gradient, Divergence, Curl, VectorZero, Laplacian
 from sympy.printing.repr import srepr
+from sympy.vector.operators import _christoffel_symbol_2nd_kind
+from sympy import zeros, Matrix, sin, cos
+
 
 R = CoordSys3D('R')
 s1 = R.x*R.y*R.z  # type: ignore
@@ -41,3 +44,41 @@ def test_Laplacian():
     assert Laplacian(v3).doit() == 2*R.i + 2*R.j + 2*R.k
     assert srepr(Laplacian(s3)) == \
             'Laplacian(Add(Pow(R.x, Integer(2)), Pow(R.y, Integer(2)), Pow(R.z, Integer(2))))'
+
+
+def test_christoffel_symbol_2nd_kind_cylindrical_system():
+    C = CoordSys3D("C", transformation="cylindrical")
+    r, theta, z = q = C.base_scalars()
+    h = C.lame_coefficients()
+
+    Gamma_r = zeros(3)
+    Gamma_theta = zeros(3)
+    Gamma_z = zeros(3)
+    for i in range(3):
+        for j in range(3):
+            Gamma_r[i, j] = _christoffel_symbol_2nd_kind(h, q, i, j, 0).doit()
+            Gamma_theta[i, j] = _christoffel_symbol_2nd_kind(h, q, i, j, 1).doit()
+            Gamma_z[i, j] = _christoffel_symbol_2nd_kind(h, q, i, j, 2).doit()
+
+    assert Gamma_r == Matrix([[0, 0, 0], [0, -r, 0], [0, 0, 0]])
+    assert Gamma_theta == Matrix([[0, 1/r, 0], [1/r, 0, 0], [0, 0, 0]])
+    assert Gamma_z == zeros(3)
+
+
+def test_christoffel_symbol_2nd_kind_spherical_system():
+    S = CoordSys3D("C", transformation="spherical")
+    r, theta, phi = q = S.base_scalars()
+    h = S.lame_coefficients()
+
+    Gamma_r = zeros(3)
+    Gamma_theta = zeros(3)
+    Gamma_phi = zeros(3)
+    for i in range(3):
+        for j in range(3):
+            Gamma_r[i, j] = _christoffel_symbol_2nd_kind(h, q, i, j, 0).doit()
+            Gamma_theta[i, j] = _christoffel_symbol_2nd_kind(h, q, i, j, 1).doit()
+            Gamma_phi[i, j] = _christoffel_symbol_2nd_kind(h, q, i, j, 2).doit()
+
+    assert Gamma_r == Matrix([[0, 0, 0], [0, -r, 0], [0, 0, -r*sin(theta)**2]])
+    assert Gamma_theta == Matrix([[0, 1/r, 0], [1/r, 0, 0], [0, 0, -sin(theta)*cos(theta)]])
+    assert Gamma_phi == Matrix([[0, 0, 1/r], [0, 0, cos(theta)/sin(theta)], [1/r, cos(theta)/sin(theta), 0]])