Fix multipatch boundary names (#880)

nutils/__init__.py

@@ -1,4 +1,4 @@
 'Numerical Utilities for Finite Element Analysis'
 
-__version__ = version = '9a33'
+__version__ = version = '9a34'
 version_name = 'jook-sing'

nutils/mesh.py

@@ -88,7 +88,7 @@ def rectilinear(richshape: Sequence[Union[int, Sequence[float]]], periodic: Sequ
 
 
 @log.withcontext
-def multipatch(patches, nelems, patchverts=None):
+def multipatch(patches, nelems, patchverts=None, space: str = 'X'):
     '''multipatch rectilinear mesh generator
 
     Generator for a :class:`~nutils.topology.MultipatchTopology` and geometry.
@@ -265,7 +265,10 @@ def multipatch(patches, nelems, patchverts=None):
                 raise ValueError('duplicate number of elements specified for patch {} in dimension {}'.format(i, dim))
             shape.append(nelems_sides[0])
         # create patch topology
-        topos.append(rectilinear(shape, root=transform.Index(ndims, i))[0])
+        topos.append(topology.StructuredTopology(space,
+            root=transform.Index(ndims, i),
+            axes=[transformseq.DimAxis(i=0, j=n, mod=0, isperiodic=False) for n in shape]))
+
         # compute patch geometry
         patchcoords = [numpy.linspace(0, 1, n+1) for n in shape]
         patchcoords = numeric.meshgrid(*patchcoords).reshape(ndims, -1)

nutils/topology.py

@@ -1997,10 +1997,8 @@ def boundary(self):
         btopos = [StructuredTopology(self.space, root=self.root, axes=self.axes[:idim] + (bndaxis,) + self.axes[idim+1:], nrefine=self.nrefine, bnames=self._bnames)
                   for idim, axis in enumerate(self.axes)
                   for bndaxis in axis.boundaries(nbounds)]
-        if not btopos:
-            return EmptyTopology(self.space, self.transforms.todims, self.ndims-1)
         bnames = [bname for bnames, axis in zip(self._bnames, self.axes) if axis.isdim and not axis.isperiodic for bname in bnames]
-        return DisjointUnionTopology(btopos, bnames)
+        return disjoint_union_topology(self.space, self.transforms.todims, self.ndims-1, btopos, bnames)
 
     @cached_property
     def interfaces(self):
@@ -2533,6 +2531,16 @@ def refined(self):
         return UnionTopology([topo.refined for topo in self._topos], self._names)
 
 
+def disjoint_union_topology(space: str, todims: int, fromdims: int, topos: Sequence[Topology], names: Sequence[str] = ()):
+    assert all(topo.space == space and topo.ndims == fromdims for topo in topos)
+    if not topos:
+        return EmptyTopology(space, todims, fromdims)
+    elif len(topos) == 1 and not names:
+        return topos[0]
+    else:
+        return DisjointUnionTopology(topos, names)
+
+
 class DisjointUnionTopology(TransformChainsTopology):
     'grouped topology'
 
@@ -2554,13 +2562,7 @@ def __init__(self, topos: Sequence[Topology], names: Sequence[str] = ()):
 
     def get_groups(self, *groups: str) -> TransformChainsTopology:
         topos = (topo if name in groups else topo.get_groups(*groups) for topo, name in itertools.zip_longest(self._topos, self._names))
-        topos = tuple(filter(None, topos))
-        if len(topos) == 0:
-            return self.empty_like()
-        elif len(topos) == 1:
-            return topos[0]
-        else:
-            return DisjointUnionTopology(topos)
+        return disjoint_union_topology(self.space, self.transforms.todims, self.ndims, tuple(filter(None, topos)))
 
     @cached_property
     def refined(self):
@@ -3055,6 +3057,7 @@ def __init__(self, topos: Sequence[Topology], connectivity):
         self._topos = tuple(topos)
         self._connectivity = numpy.array(connectivity)
 
+        assert all(isinstance(topo, StructuredTopology) for topo in self._topos)
         space = self._topos[0].space
         assert all(topo.space == space for topo in self._topos)
 
@@ -3065,10 +3068,11 @@ def __init__(self, topos: Sequence[Topology], connectivity):
             transformseq.chain([topo.opposites for topo in self._topos], self._topos[0].transforms.todims, self._topos[0].ndims))
 
         sides = {}
-        for topo, verts in zip(self._topos, self._connectivity):
+        for ipatch, verts in enumerate(self._connectivity):
             for idim, iside, idx in self._iter_boundaries():
                 bverts = verts[idx]
-                sides.setdefault(frozenset(bverts.flat), []).append((bverts, (topo, idim, iside)))
+                bname = self._topos[ipatch]._bnames[idim][iside]
+                sides.setdefault(frozenset(bverts.flat), []).append((bverts, (ipatch, bname)))
 
         self._boundaries = []
         self._interfaces = []
@@ -3079,20 +3083,14 @@ def __init__(self, topos: Sequence[Topology], connectivity):
             else:
                 if not all((bverts == bverts0).all() for bverts in other_bverts):
                     raise NotImplementedError('patch interfaces must have the same order of axes and the same orientation per axis')
-                self._interfaces.append((bverts0, data))
+                self._interfaces.append(data)
 
     def _iter_boundaries(self):
         return ((idim, iside, (slice(None),)*idim + (iside,)) for idim in range(self.ndims) for iside in (-1, 0))
 
     def get_groups(self, *groups: str) -> TransformChainsTopology:
         topos = (topo if 'patch{}'.format(i) in groups else topo.get_groups(*groups) for i, topo in enumerate(self._topos))
-        topos = tuple(filter(None, topos))
-        if len(topos) == 0:
-            return self.empty_like()
-        elif len(topos) == 1:
-            return topos[0]
-        else:
-            return DisjointUnionTopology(topos)
+        return disjoint_union_topology(self.space, self.transforms.todims, self.ndims, tuple(filter(None, topos)))
 
     def basis_std(self, degree, patchcontinuous=True):
         return self.basis_spline(degree, patchcontinuous, continuity=0)
@@ -3207,16 +3205,12 @@ def basis_patch(self):
     def boundary(self):
         'boundary'
 
-        if not self._boundaries:
-            return EmptyTopology(self.space, self.transforms.todims, self.ndims-1)
-
         subtopos = []
         subnames = []
-        for i, (topo, idim, iside) in enumerate(self._boundaries):
-            name = topo._bnames[idim][iside]
-            subtopos.append(topo.boundary[name])
-            subnames.append('patch{}-{}'.format(i, name))
-        return DisjointUnionTopology(subtopos, subnames)
+        for ipatch, bname in self._boundaries:
+            subtopos.append(self._topos[ipatch].boundary[bname])
+            subnames.append(f'patch{ipatch}-{bname}')
+        return disjoint_union_topology(self.space, self.transforms.todims, self.ndims-1, subtopos, subnames)
 
     @cached_property
     def interfaces(self):
@@ -3227,23 +3221,21 @@ def interfaces(self):
         patch via ``'intrapatch'``.
         '''
 
-        intrapatchtopo = EmptyTopology(self.space, self.transforms.todims, self.ndims-1) if not self._topos else \
-            DisjointUnionTopology([topo.interfaces for topo in self._topos])
+        intrapatchtopo = disjoint_union_topology(self.space, self.transforms.todims, self.ndims-1, [topo.interfaces for topo in self._topos])
 
         btopos = []
-        bconnectivity = []
-        for bverts, patchdata in self._interfaces:
+        for patchdata in self._interfaces:
             if len(patchdata) > 2:
                 raise ValueError('Cannot create interfaces of multipatch topologies with more than two interface connections.')
-            boundaries = [topo.boundary[topo._bnames[idim][iside]] for topo, idim, iside in patchdata]
+            boundaries = [self._topos[ipatch].boundary[bname] for ipatch, bname in patchdata]
             transforms, opposites = [btopo.transforms for btopo in boundaries]
             references, opposite_references = [btopo.references for btopo in boundaries]
             assert references == opposite_references
             # create structured topology of joined element pairs
             btopos.append(TransformChainsTopology(self.space, references, transforms, opposites))
-            bconnectivity.append(bverts)
         # create multipatch topology of interpatch boundaries
-        interpatchtopo = MultipatchTopology(btopos, bconnectivity)
+
+        interpatchtopo = disjoint_union_topology(self.space, self.transforms.todims, self.ndims-1, btopos)
 
         return DisjointUnionTopology((intrapatchtopo, interpatchtopo), ('intrapatch', 'interpatch'))
 

tests/test_topology.py

@@ -1104,10 +1104,10 @@ class multipatch_L(TestCase):
 
     def setUp(self):
         # 2---5
-        # |   |
-        # 1---4------7
-        # |   |      |
-        # 0---3------6
+        # | 1 |
+        # 1---4-------7
+        # | 0 |   2   |
+        # 0---3-------6
 
         super().setUp()
         self.domain, self.geom = mesh.multipatch(
@@ -1164,6 +1164,26 @@ def test_connectivity(self):
             self.assertEqual(trans1, interfaces2.transforms[i2])
             self.assertEqual(opp1, interfaces2.opposites[i2])
 
+    def assertCentroid(self, topo, x, y):
+        J = function.J(self.geom)
+        CA, A = topo.integrate([self.geom * J, J], degree=1)
+        cx, cy = CA / A
+        self.assertAlmostEqual(x, cx)
+        self.assertAlmostEqual(y, cy)
+
+    def test_groups(self):
+        self.assertCentroid(self.domain['patch0'], 0.5, 0.5)
+        self.assertCentroid(self.domain['patch1'], 0.5, 1.5)
+        self.assertCentroid(self.domain['patch2'], 2, 0.5)
+        self.assertCentroid(self.domain.boundary['patch0-bottom'], 0.5, 0)
+        self.assertCentroid(self.domain.boundary['patch0-left'], 0, 0.5)
+        self.assertCentroid(self.domain.boundary['patch1-left'], 0, 1.5)
+        self.assertCentroid(self.domain.boundary['patch1-top'], .5, 2)
+        self.assertCentroid(self.domain.boundary['patch1-right'], 1, 1.5)
+        self.assertCentroid(self.domain.boundary['patch2-top'], 2, 1)
+        self.assertCentroid(self.domain.boundary['patch2-right'], 3, .5)
+        self.assertCentroid(self.domain.boundary['patch2-bottom'], 2, 0)
+
 
 class multipatch_misc(TestCase):