title: How to get meshgrid of vectors of matrices in numpy?
body: <p>A normal meshgrid on two 1d vectors returns a matrix for each, containing duplicates of itself to fit the length of the other.</p>
  <pre><code>import numpy as np
  
  a, b = np.meshgrid(np.arange(2), np.arange(3, 6))
  a
  Out[22]: 
  array([[0, 1],
         [0, 1],
         [0, 1]])
  b
  Out[23]: 
  array([[3, 3],
         [4, 4],
         [5, 5]])
  
  </code></pre>
  <p>I want the same, but with each element being a n-d volume, with the meshgrid only on the 1st dimension:</p>
  <pre><code>v1
  Out[17]: 
  array([[0, 1, 2, 3, 4],
         [5, 6, 7, 8, 9]])
  
  v2
  Out[18]: 
  array([[11, 12, 13, 14, 15, 16],
         [17, 18, 19, 20, 21, 22],
         [23, 24, 25, 26, 27, 28]])
  
  v1.shape
  Out[19]: (2, 5)
  
  v2.shape
  Out[20]: (3, 6)
  
  </code></pre>
  <p>I want two meshgrids <code>v1_mesh.shape==(2, 3, 5)</code> and <code>v2_mesh.shape==(2, 3, 6)</code>.</p>
  <p><code>v1_mesh[i, :, :] == v1</code> and <code>v2_mesh[:, j, :] == v2</code> for all relevant indices, just like a standard meshgrid.</p>
  <p>That is a total of <code>2*3=6 == np.prod([v1.shape[0], v2.shape[0]])</code> combinations.</p>
  <hr />
  <p>Using <code>a, b = np.meshgrid(v1, v2)</code> gives <code>a.shape == b.shape == (np.prod(v1.shape), np.prod(v1.shape))</code> which is more combinations than I wanted. I only want combinations along the 1st axis.</p>