docs: def disvec(self, x, y, aq=None):
          """Discharge vector at `x`, `y`
          
          Returns
          -------
          
          qxqy : array size (2, naq)
              first row is Qx in each aquifer layer, second row is Qy
          """
          
          if aq is None: aq = self.aq.find_aquifer_data(x, y)
          rv = np.zeros((2, aq.naq))
          for e in aq.elementlist:
              rv += e.disvec(x, y, aq)
          rv = np.sum(rv[:, np.newaxis, :] * aq.eigvec, 2)
          return rv