diff --git a/nutils/topology.py b/nutils/topology.py index 105e8dfe5..a28dff5df 100644 --- a/nutils/topology.py +++ b/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,15 +2359,10 @@ 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' @@ -2362,6 +2370,17 @@ def __init__(self, space: str, simplices: numpy.ndarray, transforms: transformse 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) diff --git a/tests/test_topology.py b/tests/test_topology.py index 5d1308edf..ab4e0920a 100644 --- a/tests/test_topology.py +++ b/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)