Merge pull request #699 from evalf/topology

Topology

nutils/mesh.py

@@ -641,7 +641,6 @@ def unitsquare(nelems, etype):
 
     if etype == 'square':
         topo, geom = rectilinear([nelems, nelems], space=space)
-        return topo, geom / nelems
 
     elif etype in ('triangle', 'mixed'):
         simplices = numpy.concatenate([
@@ -669,15 +668,15 @@ def unitsquare(nelems, etype):
             transforms = transformseq.IndexTransforms(2, len(connectivity))
             topo = topology.ConnectedTopology(space, References.from_iter(references, 2), transforms, transforms, connectivity)
 
-        geom = (topo.basis('std', degree=1) * coords.T).sum(-1)
+        geom = topo.basis('std', degree=1) @ coords
         x, y = topo.boundary.sample('_centroid', None).eval(geom).T
-        bgroups = dict(left=x < .1, right=x > nelems-.1, bottom=y < .1, top=y > nelems-.1)
-        topo = topo.withboundary(**{name: topo.boundary[numpy.where(mask)[0]] for name, mask in bgroups.items()})
-        return topo, geom / nelems
+        btopo = topo.boundary
+        topo = topo.withboundary(left=btopo[x<.1], right=btopo[x>nelems-.1], bottom=btopo[y<.1], top=btopo[y>nelems-.1])
 
     else:
         raise Exception('invalid element type {!r}'.format(etype))
 
+    return topo, geom / nelems
 
 try:
     from math import comb as _comb  # new in Python 3.8

nutils/topology.py

@@ -276,6 +276,8 @@ def __getitem__(self, item: Any) -> 'Topology':
             return topo
         elif numeric.isintarray(item) and item.ndim == 1 or isinstance(item, Sequence) and all(isinstance(i, int) for i in item):
             return self.take(item)
+        elif numeric.isboolarray(item) and item.ndim == 1 and len(item) == len(self):
+            return self.compress(item)
         else:
             raise NotImplementedError
         if not topo:
@@ -2329,6 +2331,20 @@ def __init__(self, space, simplices: _renumber, transforms: transformseq.strictt
         references = References.uniform(element.getsimplex(transforms.fromdims), len(transforms))
         super().__init__(space, references, transforms, opposites)
 
+    def take_unchecked(self, indices):
+        space, = self.spaces
+        return SimplexTopology(space, self.simplices[indices], self.transforms[indices], self.opposites[indices])
+
+    @property
+    def boundary(self):
+        space, = self.spaces
+        ielem, iedge = (self.connectivity == -1).nonzero()
+        nd = self.ndims
+        edges = numpy.arange(nd+1).repeat(nd).reshape(nd,nd+1).T[::-1]
+        simplices = self.simplices[ielem, edges[iedge].T].T
+        transforms = self.transforms.edges(self.references)[ielem * (nd+1) + iedge]
+        return SimplexTopology(space, simplices, transforms, transforms)
+
     @property
     def connectivity(self):
         nverts = self.ndims + 1

tests/test_topology.py

@@ -104,6 +104,7 @@ def test_take(self):
         self.assertFalse(self.topo.take([]))
         for ielem in range(self.desired_nelems):
             self.assertTake(self.topo.take([ielem]), [ielem])
+            self.assertTake(self.topo[[ielem]], [ielem])
 
     def test_take_invalid_indices(self):
         with self.assertRaisesRegex(ValueError, '^expected a one-dimensional array$'):
@@ -115,6 +116,7 @@ def test_compress(self):
         self.assertFalse(self.topo.compress([False]*self.desired_nelems))
         for ielem in range(self.desired_nelems):
             self.assertTake(self.topo.compress([i == ielem for i in range(self.desired_nelems)]), [ielem])
+            self.assertTake(self.topo[numpy.arange(self.desired_nelems) == ielem], [ielem])
 
     def test_slice_invalid_dim(self):
         with self.assertRaisesRegex(IndexError, '^dimension index out of range$'):
@@ -1303,6 +1305,31 @@ def test_get_groups_parent(self):
         self.assertEqual(len(self.topo.get_groups('a')), 2)
 
 
+class SimplexTopology(TestCase, CommonTests, TransformChainsTests, ConformingTests):
+
+    def setUp(self):
+        super().setUp()
+        coords = numpy.array([[0,0],[0,1],[1,0],[1,1],[.5,.5]])
+        simplices = numpy.array([[0,1,4],[0,2,4],[1,3,4],[2,3,4]])
+        transforms = transformseq.IndexTransforms(2, len(simplices))
+        self.topo = topology.SimplexTopology('X', simplices, transforms, transforms)
+        self.geom = self.topo.basis('std', degree=1) @ coords
+        self.desired_nelems = 4
+        self.desired_spaces = 'X',
+        self.desired_space_dims = 2,
+        self.desired_ndims = 2
+        self.desired_volumes = [.25]*4
+        self.desired_references = [element.TriangleReference()]*4
+        self.desired_vertices = coords[simplices].tolist()
+
+    def test_boundary(self):
+        self.assertIsInstance(self.topo.boundary, topology.SimplexTopology)
+
+    def test_getitem(self):
+        self.assertIsInstance(self.topo[numpy.arange(4) < 2], topology.SimplexTopology)
+        self.assertIsInstance(self.topo[numpy.arange(2)], topology.SimplexTopology)
+
+
 class project(TestCase):
 
     def setUp(self):