simplify slice and remove midpoint for ndims>1

This patch simplifies the Reference object's slice method by changing the first
argument from a callable to a level array, making use of the edge_vertices
attribute to propagate this to the edges, and changes the midpoint logic for
elements of dimension >= 2 by selecting an existing vertex rather than
introducing a new midpoint.

examples/platewithhole.py

@@ -99,32 +99,32 @@ def test_spline(self):
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00033, places=5)
         with self.subTest('h1-error'):
-            self.assertAlmostEqual(err[1], .00672, places=5)
+            self.assertAlmostEqual(err[1], .00674, places=5)
         with self.subTest('constraints'):
             self.assertAlmostEqual64(cons, '''
                 eNpjaGBoYGBAxvrnGBow4X89g3NQFSjQwLAGq7i10Wus4k+NfM8fNWZgOGL89upc47WX0ozvXjAzPn1e
                 1TjnPACrACoJ''')
         with self.subTest('left-hand side'):
             self.assertAlmostEqual64(u, '''
-                eNpbZHbajIHhxzkGBhMgtgdi/XPypyRPvjFxO/PccPq5Vn2vcxr6luf+6xmcm2LMwLDQePf5c0bTzx8x
-                5D7vaTjnnIFhzbmlQPH5xhV39Y3vXlxtJHoh2EjvvLXR63MbgOIbjRdfrTXeecnUeO+Fn0Yrzj818j1/
-                FCh+xPjt1bnGay+lGd+9YGZ8+ryqcc55AK+AP/0=''')
+                eNpbb3bMjIHhxzkGBhMgtgdi/XMqp8RPvjLxOPPCcNq5Fn3Pcxr6luf+6xmcm2LMwLDQePf5c0bTzx8x
+                5DnvaTjnnIFhzbmlQPH5xgvu6hvfvbjaSPRCsJHeeWuj1+c2AMU3Gi++Wmu885Kp8d4LP41WnH9q5Hv+
+                KFD8iPHbq3ON115KM757wcz49HlV45zzAL8gQC8=''')
 
     def test_mixed(self):
         err, cons, u = main(nelems=4, etype='mixed', btype='std', degree=2, traction=.1, maxrefine=2, radius=.5, poisson=.3)
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00024, places=5)
         with self.subTest('h1-error'):
-            self.assertAlmostEqual(err[1], .00739, places=5)
+            self.assertAlmostEqual(err[1], .00740, places=5)
         with self.subTest('constraints'):
             self.assertAlmostEqual64(cons, '''
                 eNpjaGDADhlwiOEU1z8HZusbgukkg5BzRJqKFRoa1oD1HzfceA5NH9FmgKC10SuwOdONpM7DxDYa77gM
                 MueoMQPDEePzV2Hic42XXmoynnQRxvc3dryQbnz3Aoj91Mj3vJnx6fOqxjnnAQzkV94=''')
         with self.subTest('left-hand side'):
             self.assertAlmostEqual64(u, '''
-                eNoNzE8og3EcBvC3uUo5rNUOnBSK9/19n0Ic0Eo5oJBmRxcaB04kUnPgoETmT2w7LVrtMBy4auMw+35/
-                7/vaykFSFEopKTnIe/jU01PPU6FNWcQIn+Or5CBfSqCGD1uDYhi7/KbW+dma5aK65gX6Y8Po8HSzZQ7y
-                vBniHyvFV9aq17V7TK42O9kwFS9YUzxhjXIcZxLCnIzjTsfxah/BMFJotjUlZYz6xYeoPqEPKaigbKhb
-                9lOj9NGa9KgtVmqJH9UT36gcp71dEr6HaVS5GS8f46AcQ9itx739SQXdBL8dRqeTo1odox35poh2yJVh
-                apEueucsRWWPgpJFoLKPNzeHC/fU+yl48pDyMi6dCFbsBNJODNu2iawOoE4PoVdP4kH/UkZeaEDaUJQG
-                zMg/DouRUg==''')
+                eNoNzEEoQ3EcB/AXVymHtdqBkyLx3v/3LTQHtHJQKKHZ0YXMQS6sSM2BcrKMqbHTotUOw4GrthzWfr//
+                e6+nHJYUyUopKSnlHT717Vvfr0cpSWCWr/FVs1GuZdHKmb6QGMYRN9Qev1irXFUVTtAfG8agb5gtc5LX
+                zQj/WDm+s3b8bsBncosZZsNUvGEt8YI1w2lcSQRrMg9Pp/FmZ2EYOfTYmnIyR+PShLi+oA8pq5DsqxoH
+                qEvGaFdG1AErtclP6pnvVYnz/u4MVj2OZrfg53OceElE3Q482p9U0d0I2FGEnRK16SQdyjfFtEOuTFOv
+                DFGDi7QsxxSSIoIPGby7Jdy4l/5PxVeGeFu4dWLYtk+Rd5JI2SaKOoh2PYVRvYi6/qWCvNKE9KMqnViR
+                fyhZkYI=''')

nutils/element.py

@@ -231,37 +231,38 @@ def trim(self, levels, maxrefine, ndivisions):
             else self.empty if numpy.less_equal(levels, 0).all() \
             else self.with_children(cref.trim(clevels, maxrefine-1, ndivisions)
                                     for cref, clevels in zip(self.child_refs, self.child_divide(levels, maxrefine))) if maxrefine > 0 \
-            else self.slice(lambda vertices: numeric.dot(numeric.poly_eval(self._linear_bernstein, vertices), levels), ndivisions)
+            else self.slice(numeric.poly_eval(self._linear_bernstein, self.vertices) @ levels, ndivisions)
 
     @property
     def _linear_bernstein(self):
         return self.get_poly_coeffs('bernstein', degree=1)
 
-    def slice(self, levelfunc, ndivisions):
+    def slice(self, levels, ndivisions):
         # slice along levelset by recursing over dimensions
 
-        levels = levelfunc(self.vertices)
-
-        assert numeric.isint(ndivisions)
-        assert len(levels) == self.nverts
+        assert len(levels) == len(self.vertices)
         if numpy.greater_equal(levels, 0).all():
             return self
         if numpy.less_equal(levels, 0).all():
             return self.empty
         assert self.ndims >= 1
 
-        refs = [edgeref.slice(lambda vertices: levelfunc(edgetrans.apply(vertices)), ndivisions) for edgetrans, edgeref in self.edges]
+        refs = [edgeref.slice(levels[edgeverts], ndivisions) for edgeverts, edgeref in zip(self.edge_vertices, self.edge_refs)]
         if sum(ref != baseref for ref, baseref in zip(refs, self.edge_refs)) < self.ndims:
             return self
         if sum(bool(ref) for ref in refs) < self.ndims:
             return self.empty
 
-        nbins = 2**ndivisions
-
         if self.ndims == 1:
 
+            # For 1D elements a midpoint is introduced through linear
+            # interpolation of the vertex levels, followed by a binning step to
+            # remove near-vertex cuts and improve robustness for topology-wide
+            # connectivity.
+
             iedge = [i for (i,), edge in zip(self.edge_vertices, self.edge_refs) if edge]
             l0, l1 = levels[iedge]
+            nbins = 2**ndivisions
             xi = numpy.round(l0/(l0-l1) * nbins)
             if xi in (0, nbins):
                 return self.empty if xi == 0 and l1 < 0 or xi == nbins and l0 < 0 else self
@@ -270,12 +271,22 @@ def slice(self, levelfunc, ndivisions):
 
         else:
 
-            clevel = levelfunc(self.centroid[_])[0]
-            select = clevel*levels <= 0 if clevel != 0 else levels != 0
-            levels = levels[select]
-            vertices = self.vertices[select]
-            xi = numpy.round(levels/(levels-clevel) * nbins)
-            midpoint = numpy.mean(vertices + (self.centroid-vertices)*(xi/nbins)[:, _], axis=0)
+            # For higher-dimensional elements, the first vertex that is newly
+            # introduced by an edge slice is selected to serve as 'midpoint'.
+            # In case no new vertices are introduced (all edges are either
+            # fully retained or fully removed) then the first vertex is
+            # selected that occurs in only one of the edges. Either situation
+            # guarantees that the selected vertex lies on the exterior hull.
+
+            for trans, edge, emap, newedge in zip(self.edge_transforms, self.edge_refs, self.edge_vertices, refs):
+                if newedge.nverts > edge.nverts:
+                    midpoint = trans.apply(newedge.vertices[edge.nverts])
+                    break
+            else:
+                count = numpy.zeros(self.nverts, dtype=int)
+                for emap, eref in zip(self.edge_vertices, refs):
+                    count[emap[eref.simplices]] += 1
+                midpoint = self.vertices[count==1][0]
 
         return MosaicReference(self, refs, midpoint)
 

tests/test_finitecell.py

@@ -388,10 +388,10 @@ def test_topos(self):
         self.assertEqual(len(self.topoB), 2)
 
     def test_boundaries(self):
-        self.assertEqual(len(self.topoA.boundary), 11)
-        self.assertEqual(len(self.topoB.boundary), 8)
-        self.assertEqual(len(self.topoA.boundary['trimmed']), 5)
-        self.assertEqual(len(self.topoB.boundary['trimmed']), 5)
+        self.assertEqual(len(self.topoA.boundary), 9)
+        self.assertEqual(len(self.topoB.boundary), 6)
+        self.assertEqual(len(self.topoA.boundary['trimmed']), 3)
+        self.assertEqual(len(self.topoB.boundary['trimmed']), 3)
 
     def test_interfaces(self):
         self.assertEqual(len(self.topoA.interfaces), 4)

tests/test_points.py

@@ -165,17 +165,16 @@ def test_weights(self):
             self.assertLess(abs(pnt.weights.sum()-exact), 1e-15)
 
     def test_points(self):
-        self.assertEqual(self.bezier.npoints, 27)
+        self.assertEqual(self.bezier.npoints, 26)
         for x in [0., .25, .5, .75, 1.]:
             for y in [0., .25, .5, .75, 1.]:
                 if x or y:
                     self.assertIn([x, y], self.bezier.coords.tolist())
         self.assertIn([0., .125], self.bezier.coords.tolist())
-        self.assertIn([.0625, .0625], self.bezier.coords.tolist())
         self.assertIn([.125, 0.], self.bezier.coords.tolist())
 
     def test_tri(self):
-        self.assertEqual(len(self.bezier.tri), 34)
+        self.assertEqual(len(self.bezier.tri), 33)
 
     def test_hull(self):
-        self.assertEqual(len(self.bezier.hull), 18)
+        self.assertEqual(len(self.bezier.hull), 17)