def vdp(vec1, vec2):
    """
    Calculate the Vector Direct Product of two vectors.

    Parameters
    ----------

    vec1, vec2 : array-like, shape=(3, x, [y, z])
        Arrays of vector values in a 1D, 2D or 3D domain.

    Returns
    -------

    tensor : ndarray, shape=(3, 3, x, [y, z])
        Tensor field resulting from direct product of the specified vectors.

        .. math::

           \\textbf{A} = \\vec{v_1} \\vec{v_2}

    References
    ----------
    http://mathworld.wolfram.com/VectorDirectProduct.html

    """

    assert vec1.shape[0] == 3, "First argument provided is not a vector field"
    assert vec2.shape[0] == 3, "Second argument provided is not a vector field"
    assert vec1.shape == vec2.shape, """
        Shapes of vectors provided do not match: {}/{}
        """.format(vec1.shape, vec2.shape)

    tensor = vec1 * vec2.reshape(3, 1, *vec2.shape[1:])
    assert tensor.shape == (3, *vec1.shape), """
        Output field has shape {}, should be {}
        """.format(tensor.shape, (3, *vec1.shape))

    return tensor