From fcff2eaf52bb4c31c30a67e847af28d4ce30f829 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 6 May 2024 14:36:29 +0200 Subject: [PATCH] fix multipatch interfaces This patch fixes the interfaces of complex multipatch configurations such as the one now added to unittest test_topology.multipatch_misc.test_partial_ring. The patch essentially entails a large cleanup of code that went defunct after removal of edge-to-edge transformations. Additionally, the basis_std method now delegates to TrimmedTopology.basis_spline rather than TransformChainsTopology.basis_std for correct handling of edge cases. --- examples/cahnhilliard.py | 28 +++--- nutils/mesh.py | 6 +- nutils/topology.py | 196 +++++++++++++-------------------------- tests/test_topology.py | 29 ++++++ 4 files changed, 109 insertions(+), 150 deletions(-) diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py index bb1916e2b..b1d5fb83e 100644 --- a/examples/cahnhilliard.py +++ b/examples/cahnhilliard.py @@ -246,22 +246,22 @@ def test_multipatchcircle(self): args = main(epsilon=parse('5cm'), mobility=parse('1μL*s/kg'), nelems=3, etype='multipatch', degree=2, timestep=parse('1h'), endtime=parse('2h')) with self.subTest('concentration'): self.assertAlmostEqual64(args['φ'], ''' - eNoNyE9Ik3EYB3BwZohodpCSIk9NQub7/n7PA5IHu7Y6xKBbgUQK0mGUMIS6jIgmJRZitc3MEi/VxQVz - XsouDcHn+7w///B22HYU/wTBOuwUmcfPp8Yhr/JfrnOJJ3iU2/g1X+YHXKMcrdnvVLYt9NwMmlTfBn2g - ZhL71OZMl214H72yu+6KmtFdpPxFbzjcu3Rl59rmafdMe1xVS7iFdk2jKdjRcb0JHxGQhriDs/oeoRTk - SMoyiJN4JVl5I7OSdLeDqJ6y2ybtj6Ahb8NDXTr+PLpRoEWy9If2qUCz9oV/n+ZsJajil9/hcvyS77JS - iSbJcJEStGxbqWC+mHp/zc7bIdttV82cv2aS/ufYuOurLATr2qPzXjX2wx3qCZ1AHMbNBBVd8rZiDXdD - v+FM0KEJdGIEXfhtHppJ7ypdpEf2yGV0zER60ziPBc2ik/9Rltpp2tituCb6c8G5YBcHGEAU21KUKXmC - ZlzAT3knjyUvn+SefJVNWTl2TXrxHyMaybw=''') + eNoNwk1I02EcB/AGnYQCQcvDBqYYVGz7P8/v93QIyhdYK8lyDktKx8ouQ0hz6EEdjOpWdIqCXrAONrAg + AhGGBSHhxe/v+f//8yVpWIGXgk5DSAjf+HzyZsKkTZv5xdf5PN01A8aYCndwgrKq2hw0bzjEtfTdmfIK + /JBX6Z2eiN6zz+Uf+bSt76g635WoFHVO59SfyFE3KxsIqU2nyZmMOF6rjclbeSTDcmNfv3RLFEEoXEIe + s/iKT3iMJ3iJaczgPUL2s4wISQVFvEBKG92oLkSSblGCEuezvEgnKKOueWXp46tczz/pMI2qW94DkzPt + 5pDJ8yx16Sqzw/M8xt+orCvObb7IB7hAazql5laDtK4zuqSGnd3lTrcuHCid9pLumh20W7JQivkfXbFL + ckq+oDV6LHzTv+8esc0iOKOeOqlIjT9oe6SMBBE906/V35N52y9JjvMPukxLaqh03DasxN0tGZD/2JaU + /MYCpvABWZzDFaQxjl60YAgFdKMT7XiFgOwBQRG+vg==''') with self.subTest('chemical-potential'): self.assertAlmostEqual64(args['η'], ''' - eNoBggF9/lU4VTgzOIo4kDhIOCY3Lzd5Ngw4xjf1N1o3YzV70II3ijZoNuIyZcpFybM1jjXxzgY28TVH - zyLKzMjGyXfICsjfx73Hl8dyMzozgstEylLKFsqEyU3I2snuyJvHZMc7yLzHKctJy6rKqci4yGDIrcqK - yvLID8royADInMcIyJ7HI8jnx13HcMdax2fHWNPez6/JHzZENjI1V8jHx7AxHMpqx2PHpMgAyDk4SDgR - OH84gThhONA3ejY/OJw3Y84NyYY2C8x1NqU2BjfJyQTKYC9rN303ATZmNnY3TTaMNmI0e8mFyf3KKjYz - Nj80DDD+Mn3MhsqCM/bLsslNyRo3BzcuNUbLEMk7yDfI0snXyWjJu817zdTK+8hRyLPJt8jYx4/HEMiy - x5wzXTNbzP02CTdQNmLKFslCNc7MRcjSx5DJZ8hiOG04JDjJN4020jBZySIwQsvxyJDInMlnyCzIxMes - x3DHmMeAx1fHVMdfx0vHS8dsxzzHl8csxz3HSceixxfIfsLEFg==''') + eNoN0F8onWEcB/AjQ1pNoZb8u8Ki0HvO+zy/F4uJ2eTG7I9W1q602toOdS4pIXZSs+zUpJN2ywXJboTE + /H7f53lf08nalZqkJldGmXYzbj+Xn2X6StM0RKf6v0rtvKS3dI8yqFOd2W5k0i3K0MfhLPNeiuQ+NVC6 + O4NK6ec33ESKYu4/PJdJjnGzJm0jIRvIT45zcnc19SVo9PPMX3FkClOYxBg+YxuMqmt5KsOyJ9UoQ1Jm + RORCbuI2DuQAh/iNfbxCCPOSr+v0kKp1e/yQKUSH99g7olZKKt8HHE956d4G/dLFasXGKUHvqIXiulnN + +7m6TEf1pm7TnaokcMMPw2nK1eeqVtUEp86WMxuei9REup1Rm/P9xk6j32XJ1tmPZs2kTML0mj6sSYVs + mxPzw6pgwNQjKgmzZKr8wVSBncaieOq1unSb3PLgkSnFChlqp096ws3y15GNdowgiif4w3fkhYzLkjzg + AV7kM74rz/gDf+M8iV0fLvAxt8mWXAGfucYq''') if __name__ == '__main__': diff --git a/nutils/mesh.py b/nutils/mesh.py index 06bbcf824..7f04f0ff1 100644 --- a/nutils/mesh.py +++ b/nutils/mesh.py @@ -279,13 +279,9 @@ def multipatch(patches, nelems, patchverts=None): ]).T coords.append(patchcoords) - # build patch boundary data - - boundarydata = topology.MultipatchTopology.build_boundarydata(patches) - # join patch topologies, geometries - topo = topology.MultipatchTopology(tuple(map(topology.Patch, topos, patches, boundarydata))) + topo = topology.MultipatchTopology(topos, patches) funcsp = topo.basis('spline', degree=1, patchcontinuous=False) geom = (funcsp * numpy.concatenate(coords, axis=1)).sum(-1) diff --git a/nutils/topology.py b/nutils/topology.py index c893076da..ecb180875 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -3047,98 +3047,45 @@ def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs): return function.PlainBasis(hbasis_coeffs, hbasis_dofs, ndofs, self.f_index, self.f_coords) -@dataclass(eq=True, frozen=True) -class PatchBoundary: - - id: Tuple[int, ...] - dim: int - side: int - reverse: Tuple[bool, ...] - transpose: Tuple[int, ...] - - def apply_transform(self, array): - return array[tuple(slice(None, None, -1) if i else slice(None) for i in self.reverse)].transpose(self.transpose) - - -@dataclass(eq=True, frozen=True) -class Patch: - - topo: Topology - verts: types.arraydata - boundaries: Tuple[PatchBoundary, ...] - - class MultipatchTopology(TransformChainsTopology): 'multipatch topology' - @staticmethod - def build_boundarydata(connectivity): - 'build boundary data based on connectivity' - - boundarydata = [] - for patch in connectivity: - ndims = len(patch.shape) - patchboundarydata = [] - for dim, side in itertools.product(range(ndims), [-1, 0]): - # ignore vertices at opposite face - verts = numpy.array(patch) - opposite = tuple({0: -1, -1: 0}[side] if i == dim else slice(None) for i in range(ndims)) - verts[opposite] = verts.max()+1 - if len(set(verts.flat)) != 2**(ndims-1)+1: - raise NotImplementedError('Cannot compute canonical boundary if vertices are used more than once.') - # reverse axes such that lowest vertex index is at first position - reverse = tuple(map(bool, numpy.unravel_index(verts.argmin(), verts.shape))) - verts = verts[tuple(slice(None, None, -1) if i else slice(None) for i in reverse)] - # transpose such that second lowest vertex connects to lowest vertex in first dimension, third in second dimension, et cetera - k = [verts[tuple(1 if i == j else 0 for j in range(ndims))] for i in range(ndims)] - transpose = tuple(sorted(range(ndims), key=k.__getitem__)) - verts = verts.transpose(transpose) - # boundarid - boundaryid = tuple(verts[..., 0].flat) - patchboundarydata.append(PatchBoundary(boundaryid, dim, side, reverse, transpose)) - boundarydata.append(tuple(patchboundarydata)) - - return boundarydata - - def __init__(self, patches: Sequence[Patch]): - assert isinstance(patches, Sequence) and all(isinstance(patch, Patch) for patch in patches), f'patches={patches!r}' - self.patches = tuple(patches) - - space = patches[0].topo.space - assert all(patch.topo.space == space for patch in patches) - - for boundaryid, patchdata in self._patchinterfaces.items(): - if len(patchdata) == 1: - continue - transposes = set() - reverses = set() - for topo, boundary in patchdata: - assert boundary.transpose[-1] == boundary.dim - transposes.add(tuple(i-1 if i > boundary.dim else i for i in boundary.transpose[:-1])) - reverses.add(boundary.reverse[:boundary.dim]+boundary.reverse[boundary.dim+1:]) - if len(transposes) != 1 or len(reverses) != 1: - raise NotImplementedError('patch interfaces must have the same order of axes and the same orientation per axis') + def __init__(self, topos: Sequence[Topology], connectivity): + + self._topos = tuple(topos) + self._connectivity = numpy.array(connectivity) + + space = self._topos[0].space + assert all(topo.space == space for topo in self._topos) super().__init__( space, - util.sum(patch.topo.references for patch in self.patches), - transformseq.chain([patch.topo.transforms for patch in self.patches], self.patches[0].topo.transforms.todims, self.patches[0].topo.ndims), - transformseq.chain([patch.topo.opposites for patch in self.patches], self.patches[0].topo.transforms.todims, self.patches[0].topo.ndims)) + util.sum(topo.references for topo in self._topos), + transformseq.chain([topo.transforms for topo in self._topos], self._topos[0].transforms.todims, self._topos[0].ndims), + 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 idim, iside, idx in self._iter_boundaries(): + bverts = verts[idx] + sides.setdefault(frozenset(bverts.flat), []).append((bverts, (topo, idim, iside))) + + self._boundaries = [] + self._interfaces = [] + for patchdata in sides.values(): + (bverts0, *other_bverts), data = zip(*patchdata) + if not other_bverts: + self._boundaries.append(data[0]) + 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)) - @cached_property - def _patchinterfaces(self): - patchinterfaces = {} - for patch in self.patches: - for boundary in patch.boundaries: - patchinterfaces.setdefault(boundary.id, []).append((patch.topo, boundary)) - return types.frozendict({ - boundaryid: tuple(data) - for boundaryid, data in patchinterfaces.items() - if len(data) > 1 - }) + 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 = (patch.topo if 'patch{}'.format(i) in groups else patch.topo.get_groups(*groups) for i, patch in enumerate(self.patches)) + 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() @@ -3147,6 +3094,9 @@ def get_groups(self, *groups: str) -> TransformChainsTopology: else: return DisjointUnionTopology(topos) + def basis_std(self, degree, patchcontinuous=True): + return self.basis_spline(degree, patchcontinuous, continuity=0) + def basis_spline(self, degree, patchcontinuous=True, knotvalues=None, knotmultiplicities=None, *, continuity=-1): '''spline from vertices @@ -3199,8 +3149,8 @@ def basis_spline(self, degree, patchcontinuous=True, knotvalues=None, knotmultip dofmap = [] dofcount = 0 commonboundarydofs = {} - for ipatch, patch in enumerate(self.patches): - # build structured spline basis on patch `patch.topo` + for ipatch, (topo, verts) in enumerate(zip(self._topos, self._connectivity)): + # build structured spline basis on patch `topo` patchknotvalues = [] patchknotmultiplicities = [] for idim in range(self.ndims): @@ -3208,7 +3158,7 @@ def basis_spline(self, degree, patchcontinuous=True, knotvalues=None, knotmultip right = tuple(1 if j == idim else slice(None) for j in range(self.ndims)) dimknotvalues = set() dimknotmultiplicities = set() - for edge in zip(patch.verts[left].flat, patch.verts[right].flat): + for edge in zip(verts[left].flat, verts[right].flat): v = knotvalues.get(edge, knotvalues.get(None, missing)) m = knotmultiplicities.get(edge, knotmultiplicities.get(None, missing)) if v is missing: @@ -3223,17 +3173,17 @@ def basis_spline(self, degree, patchcontinuous=True, knotvalues=None, knotmultip raise 'ambiguous knot multiplicities for patch {}, dimension {}'.format(ipatch, idim) patchknotvalues.extend(dimknotvalues) patchknotmultiplicities.extend(dimknotmultiplicities) - patchcoeffs, patchdofmap, patchdofcount = patch.topo._basis_spline(degree, knotvalues=patchknotvalues, knotmultiplicities=patchknotmultiplicities, continuity=continuity) + patchcoeffs, patchdofmap, patchdofcount = topo._basis_spline(degree, knotvalues=patchknotvalues, knotmultiplicities=patchknotmultiplicities, continuity=continuity) coeffs.extend(patchcoeffs) dofmap.extend(types.frozenarray(dofs+dofcount, copy=False) for dofs in patchdofmap) if patchcontinuous: # reconstruct multidimensional dof structure dofs = dofcount + numpy.arange(numpy.prod(patchdofcount), dtype=int).reshape(patchdofcount) - for boundary in patch.boundaries: + for idim, iside, idx in self._iter_boundaries(): # get patch boundary dofs and reorder to canonical form - boundarydofs = boundary.apply_transform(dofs)[..., 0].ravel() + boundarydofs = dofs[idx].ravel() # append boundary dofs to list (in increasing order, automatic by outer loop and dof increment) - commonboundarydofs.setdefault(boundary.id, []).append(boundarydofs) + commonboundarydofs.setdefault(tuple(verts[idx].flat), []).append(boundarydofs) dofcount += numpy.prod(patchdofcount) if patchcontinuous: @@ -3248,28 +3198,25 @@ def basis_spline(self, degree, patchcontinuous=True, knotvalues=None, knotmultip def basis_patch(self): 'degree zero patchwise discontinuous basis' - transforms = transformseq.PlainTransforms(tuple((patch.topo.root,) for patch in self.patches), self.ndims, self.ndims) + transforms = transformseq.PlainTransforms(tuple((topo.root,) for topo in self._topos), self.ndims, self.ndims) index = function.transforms_index(self.space, transforms) coords = function.transforms_coords(self.space, transforms) - return function.DiscontBasis([types.frozenarray([[1]], dtype=float)]*len(self.patches), index, coords) + return function.DiscontBasis([types.frozenarray([[1]], dtype=float)]*len(self._topos), index, coords) @cached_property def boundary(self): 'boundary' + if not self._boundaries: + return EmptyTopology(self.space, self.transforms.todims, self.ndims-1) + subtopos = [] subnames = [] - for i, patch in enumerate(self.patches): - for boundary in patch.boundaries: - if boundary.id in self._patchinterfaces: - continue - name = patch.topo._bnames[boundary.dim][boundary.side] - subtopos.append(patch.topo.boundary[name]) - subnames.append('patch{}-{}'.format(i, name)) - if len(subtopos) == 0: - return EmptyTopology(self.space, self.transforms.todims, self.ndims-1) - else: - return DisjointUnionTopology(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) @cached_property def interfaces(self): @@ -3280,36 +3227,23 @@ def interfaces(self): patch via ``'intrapatch'``. ''' - intrapatchtopo = EmptyTopology(self.space, self.transforms.todims, self.ndims-1) if not self.patches else \ - DisjointUnionTopology([patch.topo.interfaces for patch in self.patches]) + intrapatchtopo = EmptyTopology(self.space, self.transforms.todims, self.ndims-1) if not self._topos else \ + DisjointUnionTopology([topo.interfaces for topo in self._topos]) btopos = [] bconnectivity = [] - for boundaryid, patchdata in self._patchinterfaces.items(): + for bverts, patchdata in self._interfaces: if len(patchdata) > 2: raise ValueError('Cannot create interfaces of multipatch topologies with more than two interface connections.') - pairs = [] - references = None - for topo, boundary in patchdata: - btopo = topo.boundary[topo._bnames[boundary.dim][boundary.side]] - if references is None: - references = numeric.asobjvector(btopo.references).reshape(btopo.shape) - references = references[tuple(_ if i == boundary.dim else slice(None) for i in range(self.ndims))] - references = boundary.apply_transform(references)[..., 0] - references = tuple(references.flat) - transforms = numeric.asobjvector(btopo.transforms).reshape(btopo.shape) - transforms = transforms[tuple(_ if i == boundary.dim else slice(None) for i in range(self.ndims))] - transforms = boundary.apply_transform(transforms)[..., 0] - pairs.append(tuple(map(transform.canonical, transforms.flat))) + boundaries = [topo.boundary[topo._bnames[idim][iside]] for topo, idim, iside 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 - references = References.from_iter(references, self.ndims-1) - transforms, opposites = pairs - transforms = transformseq.PlainTransforms(transforms, self.transforms.todims, self.ndims-1) - opposites = transformseq.PlainTransforms(opposites, self.transforms.todims, self.ndims-1) btopos.append(TransformChainsTopology(self.space, references, transforms, opposites)) - bconnectivity.append(numpy.array(boundaryid).reshape((2,)*(self.ndims-1))) + bconnectivity.append(bverts) # create multipatch topology of interpatch boundaries - interpatchtopo = MultipatchTopology(tuple(map(Patch, btopos, bconnectivity, self.build_boundarydata(bconnectivity)))) + interpatchtopo = MultipatchTopology(btopos, bconnectivity) return DisjointUnionTopology((intrapatchtopo, interpatchtopo), ('intrapatch', 'interpatch')) @@ -3317,11 +3251,11 @@ def interfaces(self): def connectivity(self): connectivity = [] patchinterfaces = {} - for patch in self.patches: # len(connectivity) represents the element offset for the current patch - ielems = numpy.arange(len(patch.topo)).reshape(patch.topo.shape) + len(connectivity) - for boundary in patch.boundaries: - patchinterfaces.setdefault(boundary.id, []).append((boundary.apply_transform(ielems)[..., 0], boundary.dim * 2 + (boundary.side == 0))) - connectivity.extend(patch.topo.connectivity + len(connectivity) * numpy.not_equal(patch.topo.connectivity, -1)) + for topo, verts in zip(self._topos, self._connectivity): # len(connectivity) represents the element offset for the current patch + ielems = numpy.arange(len(topo)).reshape(topo.shape) + len(connectivity) + for idim, iside, idx in self._iter_boundaries(): + patchinterfaces.setdefault(tuple(verts[idx].flat), []).append((ielems[idx], idim * 2 + (iside == 0))) + connectivity.extend(topo.connectivity + len(connectivity) * numpy.not_equal(topo.connectivity, -1)) connectivity = numpy.array(connectivity) for patchdata in patchinterfaces.values(): if len(patchdata) > 2: @@ -3339,6 +3273,6 @@ def connectivity(self): def refined(self): 'refine' - return MultipatchTopology(tuple(Patch(patch.topo.refined, patch.verts, patch.boundaries) for patch in self.patches)) + return MultipatchTopology([topo.refined for topo in self._topos], self._connectivity) # vim:sw=4:sts=4:et diff --git a/tests/test_topology.py b/tests/test_topology.py index 4ebcbf59e..7355cc506 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -1176,6 +1176,35 @@ def test_basis_merge_dofs(self): desired[numpy.arange(16), [0, 1, 2, 3, 4, 5, 3, 6, 2, 3, 7, 8, 3, 6, 8, 9]] = 1 numpy.testing.assert_allclose(actual, desired) + def test_ring(self): + patches = (0, 1, 3, 4), (3, 4, 5, 6), (1, 2, 4, 6) + patchverts = (0, 0), (.5, 0), (1, 0), (0, .5), (.5, .5), (0, 1), (1, 1) + # 5-------6 + # | 1 / | + # 3---4 | + # | 0 | 2 | + # 0---1---2 + domain, geom = mesh.multipatch(patches=patches, patchverts=patchverts, nelems=1) + self.assertEqual(domain.connectivity.tolist(), [[1, -1, 2, -1], [-1, 0, 2, -1], [1, -1, -1, 0]]) + self.assertEqual(len(domain.interfaces['interpatch']), 3) + self.assertEqual(len(domain.basis('std', 1)), 7) + self.assertAlmostEqual(domain.integrate(function.J(geom), degree=1), 1) + + def test_partial_ring(self): + patches = (10, 12, 6, 8, 11, 1, 7, 9), (10, 12, 6, 8, 0, 1, 2, 3), (13, 15, 10, 12, 14, 5, 11, 1), (13, 15, 10, 12, 4, 5, 0, 1) + # side 1 14--11--7 + # 5---1---9 | 2 | 0 | + # | 2 | 0 | 13--10--6 + # 15--12--8 | 3 | 1 | + # | 3 | 1 | 4---0---2 + # 5---1---3 side 2 + domain, geom = mesh.multipatch(patches=patches, nelems=1) + self.assertEqual(len(domain.interfaces['interpatch']), 5) # 5th interface via 15-12-5-1 + self.assertEqual(len(domain.boundary), 14) + self.assertEqual(len(domain.basis('std', 1)), 16) + # this test fails if MultipatchTopology.basis_std delegates to + # TransformChainsTopology.basis_std rather than MultipatchTopology.basis_spline. + class multipatch_errors(TestCase):