add Topology.partition_interfaces

Author: joostvanzwieten

nutils/topology.py

@@ -1120,6 +1120,35 @@ def interfaces_spaces_unchecked(self, __spaces: FrozenSet[str]) -> 'Topology':
 
         raise NotImplementedError
 
+    def partition_interfaces(self, part_indices: Sequence[int]) -> 'Topology':
+        '''Return the interfaces between all parts of an element partition.
+
+        Given a partition of elements, this function returns the subset of
+        interfaces that face two different parts. The orientation of the
+        interfaces is arbitrary.
+
+        Parameters
+        ----------
+        part_indices : sequence or :class:`numpy.ndarray` of :class:`int`
+            For each element the index of the part the element belongs to.
+
+        Returns
+        -------
+        :class:`Topology`
+            The interfaces between all parts of the element partition, a subset
+            of :attr:`Topology.interfaces`.
+        '''
+
+        part_indices = function.Array.cast(part_indices)
+        if part_indices.dtype not in (bool, int):
+            raise ValueError(f'expected a sequence of integer part indices but got a sequence of type {part_indices.dtype}')
+        if part_indices.shape != (len(self),):
+            raise ValueError(f'expected a sequence of {len(self)} integer part indices but got an array with shape {part_indices.shape}')
+        interfaces = self.interfaces
+        part_indices = numpy.take(part_indices, self.f_index)
+        isnt_partition_interface = interfaces.elementwise_stack(part_indices == function.opposite(part_indices))
+        return interfaces.compress(~function.eval(isnt_partition_interface))
+
     def basis_discont(self, degree: int) -> function.Basis:
         'discontinuous shape functions'
 

tests/test_topology.py

@@ -627,6 +627,11 @@ def test_elementwise_stack_mul(self):
         stack = self.topo.elementwise_stack(numpy.stack([self.topo1.f_index, self.topo2.f_index]), axis=1)
         self.assertEqual(function.eval(stack).tolist(), desired.tolist())
 
+    def test_partition_interfaces(self):
+        centers = self.topo.partition_interfaces([0, 0, 1, 2, 1, 1]).sample('gauss', 0).eval(self.geom).tolist()
+        centers.sort()
+        self.assertAllAlmostEqual(centers, [[0.5, 2.0], [1.0, 0.5], [1.0, 1.5], [1.5, 1.0]])
+
 
 class NewWithGroupAliases(TestCase, CommonTests, ConformingTests):
 
@@ -1365,6 +1370,11 @@ def setUp(self):
         self.desired_references = [element.LineReference()**2]*4
         self.desired_vertices = [[[x, y] for x in X for y in Y] for X in pairwise(range(3)) for Y in pairwise(range(3))]
 
+    def test_partition_interfaces(self):
+        centers = self.topo.partition_interfaces([0, 0, 1, 2]).sample('gauss', 0).eval(self.geom).tolist()
+        centers.sort()
+        self.assertAllAlmostEqual(centers, [[1.0, 0.5], [1.0, 1.5], [1.5, 1.0]])
+
 
 class UnionTopology(TestCase, CommonTests, TransformChainsTests):