fixes

nutils/function.py

@@ -4275,8 +4275,6 @@ def get_support(self, dof):
       raise IndexError('dof out of bounds')
     return numeric.sorted_index(self._transmap, self._parent.get_support(self._dofmap[dof]), missing='mask')
 
-<<<<<<< HEAD
-
 class ProductBasis(Basis):
 
   __slots__ = '_basis1', '_basis2'
@@ -4328,7 +4326,7 @@ class WithTransformsBasis(Basis):
   def __init__(self, parent:strictbasis, transforms:transformseq.stricttransforms, trans:types.strict[TransformChain]):
     self._parent = parent
     assert len(self._parent.transforms) == len(transforms)
-    super().__init__(ndofs=parent.ndofs, transforms=transforms, trans=trans)
+    super().__init__(ndofs=parent.ndofs, transforms=transforms, ndims=parent.ndimsdomain, trans=trans)
 
   def get_support(self, dof):
     return self._parent.get_support(dof)
@@ -4348,8 +4346,11 @@ def __init__(self, bases:types.tuple[strictbasis], trans:types.strict[TransformC
     self._bases = bases
     self._dofsplits = numpy.cumsum([0, *map(len, bases)])
     self._elemsplits = numpy.cumsum([0, *(len(basis.transforms) for basis in bases)])
-    transforms = transformseq.chain((basis.transforms for basis in bases), bases[0].transforms.fromdims)
-    super().__init__(ndofs=self._dofsplits[-1], transforms=transforms, trans=trans)
+    ndims = bases[0].ndimsdomain
+    if not all(basis.ndimsdomain == ndims for basis in bases):
+      raise ValueError
+    transforms = transformseq.chain((basis.transforms for basis in bases), bases[0].transforms.todims)
+    super().__init__(ndofs=self._dofsplits[-1], transforms=transforms, ndims=ndims, trans=trans)
 
   def get_support(self, dof):
     if numeric.isint(dof):

nutils/topology.py

@@ -406,9 +406,10 @@ def trim(self, levelset, maxrefine, ndivisions=8, name='trimmed', leveltopo=None
   @log.withcontext
   @types.apply_annotations
   def partition(self, levelset:function.asarray, maxrefine:types.strictint, posname:types.strictstr, negname:types.strictstr, *, ndivisions=8, arguments=None):
+    partsroot = function.Root('parts', 0)
     pos = self._trim(levelset, maxrefine=maxrefine, ndivisions=ndivisions, arguments=arguments)
     refs = tuple((pref, bref-pref) for bref, pref in zip(self.references, pos))
-    return PartitionedTopology(self, refs, (posname, negname))
+    return PartitionedTopology(self, partsroot, refs, (posname, negname))
 
   def subset(self, topo, newboundary=None, strict=False):
     'intersection'
@@ -1768,27 +1769,30 @@ def __init__(self, topos:types.tuple[stricttopology], names:types.tuple[types.st
       transformseq.chain((topo.transforms for topo in self._topos), tuple(root.ndims for root in roots)),
       transformseq.chain((topo.opposites for topo in self._topos), tuple(root.ndims for root in roots)))
 
-  def getitem(self, item):
-    topos = [topo if name == item else topo.getitem(item) for topo, name in itertools.zip_longest(self._topos, self._names)]
-    topos = [topo for topo in topos if not isinstance(topo, EmptyTopology)]
+  def _mk(self, topos, names=()):
     if len(topos) == 0:
       return EmptyTopology(self.roots, self.ndims)
     elif len(topos) == 1:
       return topos[0]
     else:
-      return DisjointUnionTopology(topos)
+      return DisjointUnionTopology(topos, names)
+
+  def getitem(self, item):
+    topos = [topo if name == item else topo.getitem(item) for topo, name in itertools.zip_longest(self._topos, self._names)]
+    topos = [topo for topo in topos if not isinstance(topo, EmptyTopology)]
+    return self._mk(topos)
 
   @property
   def refined(self):
-    return DisjointUnionTopology([topo.refined for topo in self._topos], self._names)
+    return self._mk([topo.refined for topo in self._topos], self._names)
 
   @property
   def boundary(self):
-    return DisjointUnionTopology([topo.boundary for topo in self._topos])
+    return self._mk([topo.boundary for topo in self._topos])
 
   @property
   def interfaces(self):
-    return DisjointUnionTopology([topo.interfaces for topo in self._topos])
+    return self._mk([topo.interfaces for topo in self._topos])
 
   def sample(self, ischeme, degree):
     transforms = self.transforms,
@@ -1800,7 +1804,7 @@ def basis(self, name, *args, **kwargs):
     if name == 'discont':
       return super().basis(name, *args, **kwargs)
     else:
-      return function.DisjointUnionBasis(topo.basis(name, *args, **kwargs) for topo in self._topos)
+      return function.DisjointUnionBasis(tuple(topo.basis(name, *args, **kwargs) for topo in self._topos), function.SelectChain(self.roots))
 
 class SubsetTopology(Topology):
   'trimmed'
@@ -2458,39 +2462,44 @@ class WithIdentifierTopology(Topology):
       :class:`nutils.transform.Identifier`.
   '''
 
+  __slots__ = '_parent', '_root', '_identifier'
+
   @types.apply_annotations
-  def __init__(self, parent:stricttopology, token):
+  def __init__(self, parent:stricttopology, root:function.strictroot, identifier:transformseq.stricttransforms):
+    assert len(identifier) == 1 and sum(identifier.todims) == 0
     self._parent = parent
-    self._token = token
-    super().__init__(parent.references,
-                     transformseq.WithIdentifierTransforms(parent.transforms, token),
-                     transformseq.WithIdentifierTransforms(parent.opposites, token))
+    self._root = root
+    self._identifier = identifier
+    super().__init__(parent.roots+(root,),
+                     parent.references,
+                     parent.transforms*identifier,
+                     parent.opposites*identifier)
 
   def basis(self, *args, **kwargs):
-    return function.WithTransformsBasis(self._parent.basis(*args, **kwargs), self.transforms)
+    return function.WithTransformsBasis(self._parent.basis(*args, **kwargs), self.transforms, function.SelectChain(self.roots))
 
   @property
   def refined(self):
-    return WithIdentifierTopology(self._parent.refined, self._token)
+    return WithIdentifierTopology(self._parent.refined, self._root, self._identifier.refined(elementseq.asreferences([element.PointReference()], 0)))
 
   @property
   def boundary(self):
-    return WithIdentifierTopology(self._parent.boundary, self._token)
+    return WithIdentifierTopology(self._parent.boundary, self._root, self._identifier)
 
   @property
   def interfaces(self):
-    return WithIdentifierTopology(self._parent.interfaces, self._token)
+    return WithIdentifierTopology(self._parent.interfaces, self._root, self._identifier)
 
   def getitem(self, item):
-    return WithIdentifierTopology(self._parent.getitem(item), self._token)
+    return WithIdentifierTopology(self._parent.getitem(item), self._root, self._identifier)
 
 class PartitionedTopology(DisjointUnionTopology):
 
-  __slots__ = 'basetopo', 'refs', 'names', 'nparts', '_parts'
+  __slots__ = 'basetopo', 'refs', 'names', 'nparts', 'partsroot', '_parts', '_partstransforms', '_nrefined'
   __cache__ = 'boundary', 'interfaces', 'refined'
 
   @types.apply_annotations
-  def __init__(self, basetopo:stricttopology, refs:types.tuple[types.tuple[element.strictreference]], names:types.tuple[types.strictstr]):
+  def __init__(self, basetopo:stricttopology, partsroot:function.strictroot, refs:types.tuple[types.tuple[element.strictreference]], names:types.tuple[types.strictstr], *, _nrefined=0):
     if len(refs) != len(basetopo):
       raise ValueError('Expected {} refs tuples but got {}.'.format(len(basetopo), len(refs)))
     self.nparts = len(refs[0]) if refs else len(names)
@@ -2502,14 +2511,19 @@ def __init__(self, basetopo:stricttopology, refs:types.tuple[types.tuple[element
       raise ValueError('Names may not contain colons.')
     if self.nparts == 0:
       raise ValueError('A partition consists of at least one part, but got zero.')
-    assert all(functools.reduce(operator.or_, prefs) == bref for bref, prefs in zip(basetopo.references, refs)), 'not a partition: union of parts is smaller then base'
+    assert all(functools.reduce(operator.or_, prefs) == bref for bref, prefs in zip(basetopo.references, refs)), 'not a partition: union of parts is smaller than base'
 
     self.basetopo = basetopo
     self.refs = refs
     self.names = names
+    self.partsroot = partsroot
+    self._nrefined = _nrefined
 
+    self._partstransforms = transformseq.IdentifierTransforms(0, partsroot.name, self.nparts)
+    for i in range(_nrefined):
+      self._partstransforms = self._partstransforms.refined(elementseq.asreferences([element.PointReference()], 0))
     indices = tuple(types.frozenarray(numpy.where(list(map(bool, prefs)))[0]) for prefs in zip(*refs))
-    self._parts = tuple(WithIdentifierTopology(SubsetTopology(basetopo, prefs), name) for name, prefs in zip(names, zip(*refs)))
+    self._parts = tuple(WithIdentifierTopology(SubsetTopology(basetopo, prefs), partsroot, self._partstransforms[i:i+1]) for i, prefs in enumerate(zip(*refs)))
     super().__init__(self._parts, names)
 
   def getitem(self, item):
@@ -2518,16 +2532,17 @@ def getitem(self, item):
     else:
       topo = self.basetopo.getitem(item)
       refs = tuple(tuple(ref & bref for ref in self.refs[self.basetopo.transforms.index(trans)]) for bref, trans in zip(topo.references, topo.transforms))
-      return PartitionedTopology(topo, refs, self.names)
+      return PartitionedTopology(topo, self.partsroot, refs, self.names)
 
   @property
   def boundary(self):
     baseboundary = self.basetopo.boundary
     brefs = []
     for bref, btrans in zip(baseboundary.references, baseboundary.transforms):
       ielem, etrans = self.basetopo.transforms.index_with_tail(btrans)
-      brefs.append(tuple(pref.get_from_trans(etrans) for pref in self.refs[ielem]))
-    return PartitionedTopology(baseboundary, brefs, self.names)
+      todims = tuple(t[-1].fromdims for t in self.basetopo.transforms[ielem])
+      brefs.append(tuple(pref.edge_refs[transform.index_edge_transforms(pref.edge_transforms, etrans, todims)] for pref in self.refs[ielem]))
+    return PartitionedTopology(baseboundary, self.partsroot, brefs, self.names, _nrefined=self._nrefined)
 
   @property
   def interfaces(self):
@@ -2537,8 +2552,9 @@ def interfaces(self):
     for ieelem, (eref, etrans, oppetrans) in enumerate(zip(baseifaces.references, baseifaces.transforms, baseifaces.opposites)):
       ielem, tail = self.basetopo.transforms.index_with_tail(etrans)
       ioppelem, opptail = self.basetopo.transforms.index_with_tail(oppetrans)
-      erefs = tuple(filter(lambda item: item[1], ((i, ref.get_from_trans(tail)) for i, ref in zip(self.names, self.refs[ielem]))))
-      opperefs = tuple(filter(lambda item: item[1], ((i, ref.get_from_trans(opptail)) for i, ref in zip(self.names, self.refs[ioppelem]))))
+      todims = tuple(t[-1].fromdims for t in self.basetopo.transforms[ielem])
+      erefs = tuple(filter(lambda item: item[1], ((i, ref.edge_refs[transform.index_edge_transforms(ref.edge_transforms, tail, todims)]) for i, ref in zip(self.names, self.refs[ielem]))))
+      opperefs = tuple(filter(lambda item: item[1], ((i, ref.edge_refs[transform.index_edge_transforms(ref.edge_transforms, opptail, todims)]) for i, ref in zip(self.names, self.refs[ioppelem]))))
       checkeref = eref.empty
       for aname, aeref in erefs:
         for bname, beref in opperefs:
@@ -2550,50 +2566,63 @@ def interfaces(self):
       assert checkeref == eref
     baseindices = {p: types.frozenarray(i, dtype=int) for p, i in baseindices.items()}
 
+    newedges = {}
+    def addnewedge(ielem, etrans):
+      edges = newedges.setdefault(ielem, [])
+      assert etrans not in edges
+      iedge = len(edges)
+      edges.append(etrans)
+      return ielem, iedge
     newreferences = {(a, b): [] for i, a in enumerate(self.names) for b in self.names[i+1:]}
     newtransforms = {(a, b): [] for i, a in enumerate(self.names) for b in self.names[i+1:]}
     newopposites = {(a, b): [] for i, a in enumerate(self.names) for b in self.names[i+1:]}
-    for baseref, partrefs, basetrans in zip(self.basetopo.references, self.refs, self.basetopo.transforms):
+    for ibase, (baseref, partrefs, basetrans) in enumerate(zip(self.basetopo.references, self.refs, self.basetopo.transforms)):
+      todims = tuple(t[-1].fromdims for t in basetrans)
       pool = {}
       for aname, aref in zip(self.names, partrefs):
         if not aref:
           continue
         for aetrans, aeref in aref.edges[baseref.nedges:]:
           if not aeref:
             continue
-          points = types.frozenarray(aetrans.apply(aeref.getpoints('bezier', 2).coords))
-          bname, beref, betrans = pool.pop((points, not aetrans.isflipped), (None, None, None))
+          points = types.frozenarray(aetrans.apply(aeref.getpoints('bezier', 2).coords), copy=False)
+          bname, beref, betrans = pool.pop(points, (None, None, None))
           if beref is None:
-            pool[(points, aetrans.isflipped)] = aname, aeref, aetrans
+            pool[points] = aname, aeref, aetrans
           else:
             assert aname != bname, 'elements are not supposed to count internal interfaces as edges'
-            assert aeref == beref
-            atrans = basetrans + (aetrans, transform.Identifier(self.ndims-1, aname))
-            btrans = basetrans + (betrans, transform.Identifier(self.ndims-1, bname))
+            # assert aeref == beref # disabled: aeref.trans is beref.trans.flipped if aeref is a ManifoldReference
             if self.names.index(aname) <= self.names.index(bname):
               iface = aname, bname
             else:
               iface = bname, aname
-              atrans, btrans = btrans, atrans
+              aetrans, betrans = betrans, aetrans
             newreferences[iface].append(aeref)
-            newtransforms[iface].append(atrans)
-            newopposites[iface].append(btrans)
+            newtransforms[iface].append(addnewedge(ibase, aetrans.separate(todims)))
+            newopposites[iface].append(addnewedge(ibase, betrans.separate(todims)))
+
       assert not pool, 'some interal edges have no opposites'
 
+    newielems, newedges = zip(*sorted(newedges.items(), key=lambda item: item[0]))
+    newoffsets = dict(zip(newielems, numpy.cumsum([0, *map(len, newedges)])))
+    newedges = transformseq.TrimmedEdgesTransforms(self.basetopo.transforms[numpy.asarray(newielems)], newedges)
     itopos = []
     inames = []
     for i, a in enumerate(self.names):
-      itopos.append(Topology(elementseq.asreferences(basereferences[a, a], self.ndims-1),
-                             transformseq.WithIdentifierTransforms(baseifaces.transforms[baseindices[a, a]], a),
-                             transformseq.WithIdentifierTransforms(baseifaces.opposites[baseindices[a, a]], a)))
+      itopos.append(Topology(self.roots,
+                             elementseq.asreferences(basereferences[a, a], self.ndims-1),
+                             baseifaces.transforms[baseindices[a, a]]*self._partstransforms[i:i+1],
+                             baseifaces.opposites[baseindices[a, a]]*self._partstransforms[i:i+1]))
       inames.append('{0}:{0}'.format(a))
-      for b in self.names[i+1:]:
-        base = Topology(elementseq.asreferences(basereferences[a, b] + basereferences[b, a], self.ndims-1),
-                        transformseq.WithIdentifierTransforms(transformseq.chain((baseifaces.transforms[baseindices[a, b]], baseifaces.opposites[baseindices[b, a]]), self.ndims-1), a),
-                        transformseq.WithIdentifierTransforms(transformseq.chain((baseifaces.opposites[baseindices[a, b]], baseifaces.transforms[baseindices[b, a]]), self.ndims-1), b))
-        new = Topology(elementseq.asreferences(newreferences[a, b], self.ndims-1),
-                       transformseq.PlainTransforms(newtransforms[a, b], self.ndims-1),
-                       transformseq.PlainTransforms(newopposites[a, b], self.ndims-1))
+      for j, b in enumerate(self.names[i+1:], i+1):
+        base = Topology(self.roots,
+                        elementseq.asreferences(basereferences[a, b] + basereferences[b, a], self.ndims-1),
+                        transformseq.chain((baseifaces.transforms[baseindices[a, b]], baseifaces.opposites[baseindices[b, a]]), self.basetopo.transforms.todims)*self._partstransforms[i:i+1],
+                        transformseq.chain((baseifaces.opposites[baseindices[a, b]], baseifaces.transforms[baseindices[b, a]]), self.basetopo.transforms.todims)*self._partstransforms[j:j+1])
+        newreferencesab = elementseq.asreferences(newreferences[a, b], self.ndims-1)
+        newtransformsab = newedges[numpy.fromiter((newoffsets[ielem]+iedge for ielem, iedge in newtransforms[a, b]), dtype=int)]
+        newoppositesab = newedges[numpy.fromiter((newoffsets[ielem]+iedge for ielem, iedge in newopposites[a, b]), dtype=int)]
+        new = Topology(self.roots, newreferencesab, newtransformsab*self._partstransforms[i:i+1], newoppositesab*self._partstransforms[j:j+1])
         itopos.append(DisjointUnionTopology((base, new)))
         inames.append('{}:{}'.format(a, b))
     return DisjointUnionTopology(itopos, inames)
@@ -2615,11 +2644,11 @@ def refined(self):
     refbasetopo = self.basetopo.refined
     refbindex = refbasetopo.transforms.index
     refinedrefs = [crefs for refs in self.refs for crefs in zip(*(ref.child_refs for ref in refs))]
-    indices = numpy.argsort([refbindex(trans+(ctrans,))
+    indices = numpy.argsort([refbindex(ctrans)
                              for trans, ref in zip(self.basetopo.transforms, self.references)
-                             for ctrans in ref.child_transforms])
+                             for ctrans in transform.child_transforms(trans, ref)])
     refinedrefs = tuple(map(refinedrefs.__getitem__, indices))
-    return PartitionedTopology(refbasetopo, refinedrefs, self.names)
+    return PartitionedTopology(refbasetopo, self.partsroot, refinedrefs, self.names, _nrefined=self._nrefined+1)
 
 class _SubsetOfPartitionedTopology(DisjointUnionTopology):
 

nutils/transform.py

@@ -635,6 +635,10 @@ def __init__(self, ndims:types.strictint, trans:stricttransformitem):
   def flipped(self):
     return Manifold(self.fromdims, self.trans.flipped)
 
+  @property
+  def isflipped(self):
+    return self.trans.isflipped
+
   def swapdown(self, other):
     if isinstance(other, (TensorChild, SimplexChild)):
       return ScaledUpdim(other, self), Identity(self.fromdims)