Skip to content

Commit

Permalink
fix empty boundary of SimplexTopology
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
joostvanzwieten committed Nov 28, 2023
1 parent 2c6c3e7 commit 41090cc
Show file tree
Hide file tree
Showing 2 changed files with 16 additions and 0 deletions.
3 changes: 3 additions & 0 deletions nutils/topology.py
Original file line number Diff line number Diff line change
Expand Up @@ -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'
Expand All @@ -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
Expand Down
13 changes: 13 additions & 0 deletions tests/test_topology.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 41090cc

Please sign in to comment.