fix weights of Sample.zip (#792)

Let sample `Z` be the zip of two samples `X` and `Y`. It is guaranteed
that the order in which `Z` traverses the elements of `X` is monotonic.
Not strict monotonic, however: if the points `p` of some element `i` of
`X` are divided over `n >= 2` elements of `Y`, then element `i` is
evaluated `n` times, each with a different subset `s_j` of points `p`
and weights `w`. If the subsets `(s_j)_j` are consecutive slices,
`max(s_j) < min(s_{j+1})`, then everything works fine. If not, for the
`j`-th evaluation of element `i` points `p[s_j]` is used together with
weights `w[t_j]` instead of `w[s_j]`, where `(t_j)_j` is a sequence of
consecutive slices such that the sizes of the subsets `s_j` and `t_j`
are equal. This patch fixes this problem and adds a test for this
situation.

Fixes: #791

nutils/sample.py

@@ -734,8 +734,10 @@ def get_evaluable_indices(self, ielem):
         return evaluable.Take(evaluable.Constant(self._indices), self._getslice(ielem))
 
     def get_evaluable_weights(self, ielem):
-        weights = self._samples[0].get_evaluable_weights(evaluable.Take(evaluable.Constant(self._ielems[0]), ielem))
-        return evaluable.Take(weights, self._getslice(ielem))
+        ielem0 = evaluable.Take(evaluable.Constant(self._ielems[0]), ielem)
+        slice0 = evaluable.Take(evaluable.Constant(self._ilocals[0]), self._getslice(ielem))
+        weights = self._samples[0].get_evaluable_weights(ielem0)
+        return evaluable._take(weights, slice0, axis=0)
 
 
 class _TakeElements(_TensorialSample):

tests/test_sample.py

@@ -300,6 +300,21 @@ def test_triplet(self):
 Zip(etype='triangle')
 Zip(etype='mixed')
 
+class ZipCornerCases(TestCase):
+
+    def test_reordered_reference_indices(self):
+        # The first sample in a zip typically maintains the order of the
+        # points. Here we test a situation where this is not true: the third
+        # point of `smpl` is located in the second element of `Y`, the other
+        # three in the first element of `Y`.
+        X, x = mesh.line([1, 5], space='X')
+        Y, y = mesh.line([0, 3, 4], space='Y')
+        smpl = X.sample('gauss', 6)
+        zipped = smpl.zip(Y.locate(y, smpl.eval(x) % 4, tol=1e-10))
+        self.assertAllAlmostEqual(zipped.eval(Y.f_index), [0, 0, 1, 0])
+        # Assert we get the correct weights (issue #791).
+        self.assertAllAlmostEqual(zipped.integrate(x * function.J(x)), 12)
+
 
 class TakeElements(TestCase, Common):