From c034cb93546e12dc3761a0a325f50ba26f628766 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 4 Mar 2024 10:42:04 +0100 Subject: [PATCH 1/7] fix typos in internal SI function names This patch replaces the _drop suffix to _op as it was originally intended. --- nutils/SI.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nutils/SI.py b/nutils/SI.py index aa7857865..473c0e45c 100644 --- a/nutils/SI.py +++ b/nutils/SI.py @@ -371,13 +371,13 @@ def __pow_like(op, *args, **kwargs): return (dim0**args[1]).wrap(op(arg0, *args[1:], **kwargs)) @register('isfinite', 'isnan', 'shape', 'ndim', 'size', 'normal', 'normalized') - def __unary_drop(op, *args, **kwargs): + def __unary_op(op, *args, **kwargs): (_dim0, arg0), = Quantity.__unpack(args[0]) return op(arg0, *args[1:], **kwargs) @register('lt', 'le', 'eq', 'ne', 'gt', 'ge', 'equal', 'not_equal', 'less', 'less_equal', 'greater', 'greater_equal') - def __binary_drop(op, *args, **kwargs): + def __binary_op(op, *args, **kwargs): (dim0, arg0), (dim1, arg1) = Quantity.__unpack(args[0], args[1]) if dim0 != dim1: raise TypeError(f'incompatible arguments for {op.__name__}: {dim0.__name__}, {dim1.__name__}') From fc25ae2e4df7dbfda573e959240e0d81fc708d2e Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 4 Mar 2024 11:34:03 +0100 Subject: [PATCH 2/7] fix logging of parameters without type annotation This patch moves the type interence from util.cli to the _infer_type helper function, and uses this in log_arguments. This fixes the situation in which arguments without type annotation were not logged correctly due to failing stringly serialization. --- nutils/_util.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/nutils/_util.py b/nutils/_util.py index 8fb4884f0..42471435d 100644 --- a/nutils/_util.py +++ b/nutils/_util.py @@ -519,7 +519,7 @@ def log_arguments(*args, **kwargs): with treelog.context('arguments'): for k, v in bound.arguments.items(): try: - s = stringly.dumps(sig.parameters[k].annotation, v) + s = stringly.dumps(_infer_type(sig.parameters[k]), v) except: s = str(v) treelog.info(f'{k}={s}') @@ -694,6 +694,16 @@ def add_htmllog(outrootdir: str = '~/public_html', outrooturi: str = '', scriptn treelog.info(f'log written to: {loguri}') +def _infer_type(param): + '''Infer argument type from annotation or default value.''' + + if param.annotation is not param.empty: + return param.annotation + if param.default is not param.empty: + return type(param.default) + raise Exception(f'cannot determine type for argument {param.name!r}') + + def cli(f, *, argv=None): '''Call a function using command line arguments.''' @@ -706,11 +716,7 @@ def cli(f, *, argv=None): mandatory = set() for param in inspect.signature(f).parameters.values(): - T = param.annotation - if T == param.empty and param.default != param.empty: - T = type(param.default) - if T == param.empty: - raise Exception(f'cannot determine type for argument {param.name!r}') + T = _infer_type(param) try: s = stringly.serializer.get(T) except Exception as e: From af52b214cb5f2267af86c90619677a3bd68b063f Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Wed, 13 Mar 2024 16:55:00 +0100 Subject: [PATCH 3/7] fix missing super()._simplified() in Mod --- nutils/evaluable.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index ee62106eb..d1c6b475b 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -2417,6 +2417,7 @@ def _simplified(self): lower_dividend, upper_dividend = dividend._intbounds if 0 <= lower_dividend and upper_dividend < lower_divisor: return dividend + return super()._simplified() class ArcTan2(Pointwise): From 80512850475d4bc28edbcbff035dc7dcad8c6c22 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 10 Apr 2024 15:17:22 +0200 Subject: [PATCH 4/7] fix signature handling in Evaluable._iter_stack This patch improves the parsing in _iter_stack of a function signature that contains parameters of kind VAR_POSITIONAL, by limiting itself to parameters of kind POSITIONAL_OR_KEYWORD. --- nutils/evaluable.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index d1c6b475b..7d036598b 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -415,12 +415,16 @@ def _iter_stack(self): for i, (op, indices) in enumerate(self.serialized, start=1): s = [f'%{i} = {op}'] if indices: + args = [f'%{i}' for i in indices] try: sig = inspect.signature(op.evalf) except ValueError: - s.extend(f'%{i}' for i in indices) + pass else: - s.extend(f'{param}=%{i}' for param, i in zip(sig.parameters, indices)) + for i, param in enumerate(sig.parameters.values()): + if i < len(args) and param.kind == param.POSITIONAL_OR_KEYWORD: + args[i] = param.name + '=' + args[i] + s.extend(args) yield ' '.join(s) def _format_stack(self, values, e): From 91d9e417583aa352e6c0397a93c72a33ac33de5e Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 31 Jan 2024 20:49:48 +0100 Subject: [PATCH 5/7] fix simplification of Equals involving Range This patch adds a condition to Equals.simplify to avoid that r[:,_] == r[:,_] is mistaken for r[:,_] == r[_,:], with r a Range and _ a newaxis. Though the former case of trivial equality is handled first, the bug could be triggered in case both sides of the equation are constructed in different order and therefore not identical, causing it to fall through to the second test which failed to check the insertion axes. --- nutils/evaluable.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index 7d036598b..9bd60f538 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -2452,7 +2452,7 @@ def _simplified(self): if self.ndim == 2: u1, w1 = unalign(a1) u2, w2 = unalign(a2) - if u1 == u2 and isinstance(u1, Range): + if u1.ndim == u2.ndim == 1 and u1 == u2 and w1 != w2 and isinstance(u1, Range): # NOTE: Once we introduce isunique we can relax the Range bound return Diagonalize(ones(u1.shape, bool)) return super()._simplified() From 9372fbcab6d8aab6de4dd430ca14854f416960bd Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 6 May 2024 14:36:29 +0200 Subject: [PATCH 6/7] 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 ca84cfe71..fc28219a3 100644 --- a/nutils/mesh.py +++ b/nutils/mesh.py @@ -281,13 +281,9 @@ def multipatch(patches, nelems, patchverts=None, name: Optional[str] = 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 2ea144a8e..42287f581 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -3039,98 +3039,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() @@ -3139,6 +3086,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 @@ -3191,8 +3141,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): @@ -3200,7 +3150,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: @@ -3215,17 +3165,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: @@ -3240,28 +3190,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): @@ -3272,36 +3219,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')) @@ -3309,11 +3243,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: @@ -3331,6 +3265,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 7f5a4b8bd..7f275581d 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -1170,6 +1170,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): From 1bae740d17b116c2d49e338ff4e39ce1c6ea0598 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 18 Dec 2023 12:12:15 +0100 Subject: [PATCH 7/7] Revert "set upper limit to mkl version" This reverts commit 258afac5655e98c05d9b3709efa6324dc977cde6. --- .github/workflows/test.yaml | 2 +- pyproject.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index f3f403e13..f871394db 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -102,7 +102,7 @@ jobs: - name: Configure MKL if: ${{ matrix.matrix-backend == 'mkl' }} run: | - python -um pip install --upgrade --upgrade-strategy eager 'mkl<2024' + python -um pip install --upgrade --upgrade-strategy eager mkl python -um devtools.gha.configure_mkl - name: Test run: python -um coverage run -m unittest discover -b -q -t . -s tests diff --git a/pyproject.toml b/pyproject.toml index 90d66d321..1bbd1861b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ classifiers = [ [project.optional-dependencies] docs = ["Sphinx >=1.8,<8"] export_mpl = ["matplotlib >=3.3,<4"] -matrix_mkl = ["mkl <2024"] +matrix_mkl = ["mkl"] matrix_scipy = ["scipy >=0.13,<2"] import_gmsh = ["meshio >=4,<6"]