fix simplex boundary

nutils/topology.py

@@ -1862,6 +1862,19 @@ def __or__(self, other):
     def __rsub__(self, other):
         return other
 
+    @property
+    def connectivity(self):
+        return ()
+
+    @property
+    def boundary(self):
+        if self.ndims:
+            return EmptyTopology(self.space, self.transforms.todims, self.ndims - 1)
+        else:
+            raise ValueError('A 0D topology has no boundary.')
+
+    interfaces = boundary
+
 
 def StructuredLine(space, root: transform.TransformItem, i: int, j: int, periodic: bool = False, bnames: Optional[Tuple[str, str]] = None):
     assert isinstance(i, int), f'i={i!r}'
@@ -2346,22 +2359,28 @@ def __init__(self, space: str, references: References, transforms: transformseq.
 class SimplexTopology(TransformChainsTopology):
     'simpex topology'
 
-    def _renumber(simplices):
-        simplices = numpy.asarray(simplices)
-        keep = numpy.zeros(simplices.max()+1, dtype=bool)
-        keep[simplices.flat] = True
-        return types.arraydata(simplices if keep.all() else (numpy.cumsum(keep)-1)[simplices])
-
     def __init__(self, space: str, simplices: numpy.ndarray, transforms: transformseq.Transforms, opposites: transformseq.Transforms):
         assert isinstance(space, str), f'space={space!r}'
         assert isinstance(simplices, numpy.ndarray), f'simplices={simplices!r}'
+        assert len(simplices), f'simplices is empty'
         assert simplices.shape == (len(transforms), transforms.fromdims+1)
         self.simplices = numpy.asarray(simplices)
         assert numpy.greater(self.simplices[:, 1:], self.simplices[:, :-1]).all(), 'nodes should be sorted'
         assert not numpy.equal(self.simplices[:, 1:], self.simplices[:, :-1]).all(), 'duplicate nodes'
         references = References.uniform(element.getsimplex(transforms.fromdims), len(transforms))
         super().__init__(space, references, transforms, opposites)
 
+    @cached_property
+    def contiguous_simplices(self):
+        simplices = numpy.asarray(self.simplices)
+        keep = numpy.zeros(simplices.max()+1, dtype=bool)
+        keep[simplices.flat] = True
+        return simplices if keep.all() else (numpy.cumsum(keep)-1)[simplices]
+
+    @cached_property
+    def nverts(self):
+        return self.contiguous_simplices.max() + 1
+
     def take_unchecked(self, indices):
         space, = self.spaces
         return SimplexTopology(space, self.simplices[indices], self.transforms[indices], self.opposites[indices])
@@ -2370,6 +2389,8 @@ def take_unchecked(self, indices):
     def boundary(self):
         space, = self.spaces
         ielem, iedge = (self.connectivity == -1).nonzero()
+        if len(ielem) == 0:
+            return self.empty_like().boundary
         nd = self.ndims
         edges = numpy.arange(nd+1).repeat(nd).reshape(nd,nd+1).T[::-1]
         simplices = self.simplices[ielem, edges[iedge].T].T
@@ -2394,7 +2415,7 @@ def connectivity(self):
     def basis_std(self, degree):
         if degree == 1:
             coeffs = element.getsimplex(self.ndims).get_poly_coeffs('bernstein', degree=1)
-            return function.PlainBasis([coeffs] * len(self), self.simplices, self.simplices.max()+1, self.f_index, self.f_coords)
+            return function.PlainBasis([coeffs] * len(self), self.contiguous_simplices, self.nverts, self.f_index, self.f_coords)
         return super().basis_std(degree)
 
     def basis_bubble(self):
@@ -2406,9 +2427,8 @@ def basis_bubble(self):
         coeffs[:-1] = poly.change_degree(bernstein, self.ndims, 1 + self.ndims) - bubble[None] / (self.ndims+1)
         coeffs[-1] = bubble
         coeffs = types.frozenarray(coeffs, copy=False)
-        nverts = self.simplices.max() + 1
-        ndofs = nverts + len(self)
-        nmap = [types.frozenarray(numpy.hstack([idofs, nverts+ielem]), copy=False) for ielem, idofs in enumerate(self.simplices)]
+        ndofs = self.nverts + len(self)
+        nmap = [types.frozenarray(numpy.hstack([idofs, self.nverts+ielem]), copy=False) for ielem, idofs in enumerate(self.contiguous_simplices)]
         return function.PlainBasis([coeffs] * len(self), nmap, ndofs, self.f_index, self.f_coords)
 
 

tests/test_topology.py

@@ -1352,9 +1352,30 @@ def setUp(self):
         self.desired_references = [element.TriangleReference()]*4
         self.desired_vertices = coords[simplices].tolist()
 
+    def test_contiguous_simplices(self):
+        # Same structure as the simplices in `setUp`, but renumbered.
+        simplices = numpy.array([[10,13,25],[10,18,25],[13,21,25],[18,21,25]])
+        transforms = transformseq.IndexTransforms(2, len(simplices))
+        topo = topology.SimplexTopology('X', simplices, transforms, transforms)
+        self.assertEqual(self.topo.contiguous_simplices.tolist(), self.topo.simplices.tolist())
+        self.assertEqual(self.topo.nverts, 5)
+
     def test_boundary(self):
         self.assertIsInstance(self.topo.boundary, topology.SimplexTopology)
 
+    def test_empty_boundary(self):
+        simplices = numpy.array([[0, 1, 2, 3]])
+        transforms = transformseq.IndexTransforms(3, len(simplices))
+        topo = topology.SimplexTopology('X', simplices, transforms, transforms)
+        self.assertEqual(len(topo.boundary), 4)
+        self.assertEqual(len(topo.boundary.boundary), 0)
+
+    def test_empty_interfaces(self):
+        simplices = numpy.array([[0, 1, 2, 3]])
+        transforms = transformseq.IndexTransforms(3, len(simplices))
+        topo = topology.SimplexTopology('X', simplices, transforms, transforms)
+        self.assertEqual(len(topo.interfaces), 0)
+
     def test_getitem(self):
         self.assertIsInstance(self.topo[numpy.arange(4) < 2], topology.SimplexTopology)
         self.assertIsInstance(self.topo[numpy.arange(2)], topology.SimplexTopology)