From 7f323a87fbb7761b695aac29e3d991c065ef2e57 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 12 Jun 2024 13:55:19 +0200 Subject: [PATCH 1/5] introduce disjoint_union_topology helper function This patch introduces the topology.disjoint_union_topology helper function to deduplicate logic found in several places of the topology module, to return either an EmptyTopology if the list of topologies is empty, or the topology itself if it is the only item in the list. --- nutils/topology.py | 38 +++++++++++++++----------------------- 1 file changed, 15 insertions(+), 23 deletions(-) diff --git a/nutils/topology.py b/nutils/topology.py index 42287f581..e9cff4e74 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -1989,10 +1989,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): @@ -2525,6 +2523,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' @@ -2546,13 +2554,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): @@ -3078,13 +3080,7 @@ def _iter_boundaries(self): 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) @@ -3199,16 +3195,13 @@ 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) + return disjoint_union_topology(self.space, self.transforms.todims, self.ndims-1, subtopos, subnames) @cached_property def interfaces(self): @@ -3219,8 +3212,7 @@ 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 = [] From e7c4ec252f57504f2c6908735d292650700e310f Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 12 Jun 2024 13:55:57 +0200 Subject: [PATCH 2/5] return multipatch interfaces as disjoint union This patch changes the return type of MultipatchTopology.inferfaces, from a MultipatchTopology to a DisjointUnionTopology. The reason is that the former structure added no benefits and broke the assumption that patch topologies are structured, which is now added as an assertion. --- nutils/topology.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/nutils/topology.py b/nutils/topology.py index e9cff4e74..bbd4711d6 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -3073,7 +3073,7 @@ 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)) @@ -3215,8 +3215,7 @@ def interfaces(self): 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] @@ -3225,9 +3224,9 @@ def interfaces(self): 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')) From 96cb656c5f947701aca3794ce3a7868f12fa8654 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 12 Jun 2024 15:01:23 +0200 Subject: [PATCH 3/5] limit multipatch topologies to structured patches This patch changes mesh.multipatch so that patches are always of type StructuredTopology, and adds an assertion for this in the constructor. --- nutils/mesh.py | 7 +++++-- nutils/topology.py | 1 + 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/nutils/mesh.py b/nutils/mesh.py index fc28219a3..ee2f0467a 100644 --- a/nutils/mesh.py +++ b/nutils/mesh.py @@ -87,7 +87,7 @@ def rectilinear(richshape: Sequence[Union[int, Sequence[float]]], periodic: Sequ @log.withcontext -def multipatch(patches, nelems, patchverts=None, name: Optional[str] = None): +def multipatch(patches, nelems, patchverts=None, name: Optional[str] = None, space: str = 'X'): '''multipatch rectilinear mesh generator Generator for a :class:`~nutils.topology.MultipatchTopology` and geometry. @@ -267,7 +267,10 @@ def multipatch(patches, nelems, patchverts=None, name: Optional[str] = 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) diff --git a/nutils/topology.py b/nutils/topology.py index bbd4711d6..6dde72fa2 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -3049,6 +3049,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) From 18410d922dcbd477ea90cae1293857f526492454 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 11 Jun 2024 16:59:50 +0200 Subject: [PATCH 4/5] fix group names for multipatch boundary, interface This patch fixes the "patch{i}-{bname}" names of multipatch boundary groups, in which the indices i were broken by commit fcff2eaf. --- nutils/topology.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/nutils/topology.py b/nutils/topology.py index 6dde72fa2..33a03e6cf 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -3060,10 +3060,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 = [] @@ -3198,10 +3199,9 @@ def boundary(self): 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)) + 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 @@ -3219,7 +3219,7 @@ def interfaces(self): 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 From 8cff74dbb403aff7a5c41eb8203b3e11f9e3cf26 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 12 Jun 2024 13:17:25 +0200 Subject: [PATCH 5/5] add unittest multipatch_L.test_groups This patch adds the unittest tests.test_topology.multipatch_L.test_groups to confirm that multipatch boundary group names are formed correctly. --- tests/test_topology.py | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/tests/test_topology.py b/tests/test_topology.py index 7f275581d..9185c7e27 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -1098,10 +1098,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( @@ -1158,6 +1158,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):