From 9b9fb80a245a45c58f652dc027809c9c1c7302e7 Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Tue, 28 Nov 2023 15:35:42 +0100 Subject: [PATCH 1/4] fix basis of SimplexTopology The implementations of the degree one std and the bubble bases for the `SimplexTopology` incorrectly assume that the `simplices` array contains vertex in `range(nverts)` without holes, and use this array as is as dofs, thereby introducing zeros in the basis for every dof in `set(range(simplices.max() + 1)) - set(simplices.flat)`. The boundary of a `SimplexTopology`, however, simply returns a subset of the volume `simplices` array, without renumbering. This bug was introduced by 7e29a283, which removes `_renumber` as preprocessor for the `simplices` parameter of the `SimplexTopology` constructor. To fix this problem this patch adds a `contiguous_simplices` property and uses thses as dofs. The alternative, requiring constructing a `SimplexTopology` with contiguous `simplices`, was rejected because this puts a computational burden on every boundary of a `SimplexTopology`, even if we don't create a basis. --- nutils/topology.py | 24 ++++++++++++++---------- tests/test_topology.py | 8 ++++++++ 2 files changed, 22 insertions(+), 10 deletions(-) diff --git a/nutils/topology.py b/nutils/topology.py index 105e8dfe5..54ab17e33 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -2346,12 +2346,6 @@ 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}' @@ -2362,6 +2356,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]) @@ -2394,7 +2399,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 +2411,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..dd8eb6644 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -1352,6 +1352,14 @@ 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) From 43c857f5fbfa5f230010569eefcae88f2d0e9b40 Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Tue, 28 Nov 2023 15:57:47 +0100 Subject: [PATCH 2/4] add connectivity to EmptyTopology This patch adds the `connectivity` attribute to `EmptyTopology` for completeness. --- nutils/topology.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/nutils/topology.py b/nutils/topology.py index 54ab17e33..72d377397 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -1862,6 +1862,10 @@ def __or__(self, other): def __rsub__(self, other): return other + @property + def connectivity(self): + return () + 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}' From 2c6c3e72b73039a7955542762212e7e525dabd7c Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Tue, 28 Nov 2023 16:12:29 +0100 Subject: [PATCH 3/4] add EmptyTopology.boundary,interfaces This patch implements the `boundary` and `interfaces` attributes for the `EmptyTopology`, which is needed by a subsequent patch fixing empty boundaries of a `SimplexTopology`. --- nutils/topology.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/nutils/topology.py b/nutils/topology.py index 72d377397..69c7f151b 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -1866,6 +1866,15 @@ def __rsub__(self, other): 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}' From 41090cc60be679e3d0b07a6c8b93a885bf73a878 Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Tue, 28 Nov 2023 16:15:37 +0100 Subject: [PATCH 4/4] fix empty boundary of SimplexTopology The constructor of `SimplexTopology` trips over assert numpy.greater(self.simplices[:, 1:], self.simplices[:, :-1]).all() if there are no simplices. This is triggered when taking the boundary of a fully periodic mesh. This patch fixes the problem by returning an `EmptyTopology` when the boundary is empty. In addition this patch adds a slightly improved error message when trying to create an empty `SimplexTopology. --- nutils/topology.py | 3 +++ tests/test_topology.py | 13 +++++++++++++ 2 files changed, 16 insertions(+) diff --git a/nutils/topology.py b/nutils/topology.py index 69c7f151b..a28dff5df 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -2362,6 +2362,7 @@ class SimplexTopology(TransformChainsTopology): 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' @@ -2388,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 diff --git a/tests/test_topology.py b/tests/test_topology.py index dd8eb6644..ab4e0920a 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -1363,6 +1363,19 @@ def test_contiguous_simplices(self): 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)