From 305ce6ab9779ee3347abda191a2952d339fd0c64 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 4 Jul 2022 17:24:40 +0200
Subject: [PATCH 01/19] rename nvertices_per_level to _nlinear_per_level

This patch renames the Reference object's nvertices_per_level attribute to
_nlinear_per_level to more accurately reflect its meaning.
---
 nutils/element.py  | 34 +++++++++++++++++-----------------
 nutils/topology.py |  2 +-
 2 files changed, 18 insertions(+), 18 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index d428b18d9..11333d2d5 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -169,7 +169,7 @@ def centroid(self):
     def trim(self, levels, maxrefine, ndivisions):
         'trim element along levelset'
 
-        assert len(levels) == self.nvertices_by_level(maxrefine)
+        assert len(levels) == self._nlinear_by_level(maxrefine)
         return self if not self or numpy.greater_equal(levels, 0).all() \
             else self.empty if numpy.less_equal(levels, 0).all() \
             else self.with_children(cref.trim(clevels, maxrefine-1, ndivisions)
@@ -249,7 +249,7 @@ def check_edges(self, tol=1e-15, print=print):
     def vertex_cover(self, ctransforms, maxrefine):
         if maxrefine < 0:
             raise Exception('maxrefine is too low')
-        npoints = self.nvertices_by_level(maxrefine)
+        npoints = self._nlinear_by_level(maxrefine)
         allindices = numpy.arange(npoints)
         if len(ctransforms) == 1:
             ctrans, = ctransforms
@@ -440,7 +440,7 @@ def __init__(self):
     def getpoints(self, ischeme, degree):
         return points.CoordsWeightsPoints(numpy.empty([1, 0]), [1.])
 
-    def nvertices_by_level(self, n):
+    def _nlinear_by_level(self, n):
         return 1
 
     def child_divide(self, vals, n):
@@ -461,12 +461,12 @@ def getpoints(self, ischeme, degree):
             return points.CoordsUniformPoints(numpy.arange(.5, degree)[:, _] / degree, 1)
         return super().getpoints(ischeme, degree)
 
-    def nvertices_by_level(self, n):
+    def _nlinear_by_level(self, n):
         return 2**n + 1
 
     def child_divide(self, vals, n):
         assert n > 0
-        assert len(vals) == self.nvertices_by_level(n)
+        assert len(vals) == self._nlinear_by_level(n)
         m = (len(vals)+1) // 2
         return vals[:m], vals[m-1:]
 
@@ -491,12 +491,12 @@ def getpoints(self, ischeme, degree):
             return points.CoordsUniformPoints(coords.T, .5)
         return super().getpoints(ischeme, degree)
 
-    def nvertices_by_level(self, n):
+    def _nlinear_by_level(self, n):
         m = 2**n + 1
         return ((m+1)*m) // 2
 
     def child_divide(self, vals, n):
-        assert len(vals) == self.nvertices_by_level(n)
+        assert len(vals) == self._nlinear_by_level(n)
         np = 1 + 2**n  # points along parent edge
         mp = 1 + 2**(n-1)  # points along child edge
         cvals = []
@@ -543,12 +543,12 @@ def getindices_vertex(self, n):
         indis = numpy.arange(m)
         return numpy.array([[i, j, k] for k in indis for j in indis[:m-k] for i in indis[:m-j-k]])
 
-    def nvertices_by_level(self, n):
+    def _nlinear_by_level(self, n):
         m = 2**n+1
         return ((m+2)*(m+1)*m)//6
 
     def child_divide(self, vals, n):
-        assert len(vals) == self.nvertices_by_level(n)
+        assert len(vals) == self._nlinear_by_level(n)
 
         child_indices = self.getindices_vertex(1)
 
@@ -598,12 +598,12 @@ def vertices(self):
     def centroid(self):
         return types.frozenarray(numpy.concatenate([self.ref1.centroid, self.ref2.centroid]), copy=False)
 
-    def nvertices_by_level(self, n):
-        return self.ref1.nvertices_by_level(n) * self.ref2.nvertices_by_level(n)
+    def _nlinear_by_level(self, n):
+        return self.ref1._nlinear_by_level(n) * self.ref2._nlinear_by_level(n)
 
     def child_divide(self, vals, n):
-        np1 = self.ref1.nvertices_by_level(n)
-        np2 = self.ref2.nvertices_by_level(n)
+        np1 = self.ref1._nlinear_by_level(n)
+        np2 = self.ref2._nlinear_by_level(n)
         return [v2.swapaxes(0, 1).reshape((-1,)+vals.shape[1:])
                 for v1 in self.ref1.child_divide(vals.reshape((np1, np2)+vals.shape[1:]), n)
                 for v2 in self.ref2.child_divide(v1.swapaxes(0, 1), n)]
@@ -870,8 +870,8 @@ def check_edges(self, tol=1e-15, print=print):
     def vertices(self):
         return self.baseref.vertices
 
-    def nvertices_by_level(self, n):
-        return self.baseref.nvertices_by_level(n)
+    def _nlinear_by_level(self, n):
+        return self.baseref._nlinear_by_level(n)
 
     def child_divide(self, vals, n):
         return self.baseref.child_divide(vals, n)
@@ -1062,8 +1062,8 @@ def __rsub__(self, other):
             return MosaicReference(other, inv_edge_refs, self._midpoint)
         return NotImplemented
 
-    def nvertices_by_level(self, n):
-        return self.baseref.nvertices_by_level(n)
+    def _nlinear_by_level(self, n):
+        return self.baseref._nlinear_by_level(n)
 
     @property
     def subrefs(self):
diff --git a/nutils/topology.py b/nutils/topology.py
index 2dcbae5c9..d7998e3fd 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -1438,7 +1438,7 @@ def trim(self, levelset, maxrefine, ndivisions=8, name='trimmed', leveltopo=None
             fcache = cache.WrapperCache()
             with log.iter.percentage('trimming', self.references, self.transforms, bins) as items:
                 for ref, trans, ctransforms in items:
-                    levels = numpy.empty(ref.nvertices_by_level(maxrefine))
+                    levels = numpy.empty(ref._nlinear_by_level(maxrefine))
                     cover = list(fcache[ref.vertex_cover](frozenset(ctransforms), maxrefine))
                     # confirm cover and greedily optimize order
                     mask = numpy.ones(len(levels), dtype=bool)

From f24ca96d0d3b43ac4918c5d5be40f0c08f4c05b6 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 4 Jul 2022 17:27:19 +0200
Subject: [PATCH 02/19] rename vertex_cover to _linear_cover

This patch renames the Reference object's vertex_cover property to
_linear_cover to more accurately reflect its meaning.
---
 nutils/element.py  | 4 ++--
 nutils/topology.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 11333d2d5..b619ff81b 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -246,7 +246,7 @@ def check_edges(self, tol=1e-15, print=print):
         if numpy.greater(abs(volume - self.volume), tol).any():
             print('divergence check failed: {} != {}'.format(volume, self.volume))
 
-    def vertex_cover(self, ctransforms, maxrefine):
+    def _linear_cover(self, ctransforms, maxrefine):
         if maxrefine < 0:
             raise Exception('maxrefine is too low')
         npoints = self._nlinear_by_level(maxrefine)
@@ -266,7 +266,7 @@ def vertex_cover(self, ctransforms, maxrefine):
         fcache = cache.WrapperCache()
         return tuple(((ctrans,) + trans, points, cindices[indices])
                      for ctrans, cref, cbin, cindices in zip(self.child_transforms, self.child_refs, cbins, self.child_divide(allindices, maxrefine))
-                     for trans, points, indices in fcache[cref.vertex_cover](frozenset(cbin), maxrefine-1))
+                     for trans, points, indices in fcache[cref._linear_cover](frozenset(cbin), maxrefine-1))
 
     def __str__(self):
         return self.__class__.__name__
diff --git a/nutils/topology.py b/nutils/topology.py
index d7998e3fd..eb8ffd76a 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -1439,7 +1439,7 @@ def trim(self, levelset, maxrefine, ndivisions=8, name='trimmed', leveltopo=None
             with log.iter.percentage('trimming', self.references, self.transforms, bins) as items:
                 for ref, trans, ctransforms in items:
                     levels = numpy.empty(ref._nlinear_by_level(maxrefine))
-                    cover = list(fcache[ref.vertex_cover](frozenset(ctransforms), maxrefine))
+                    cover = list(fcache[ref._linear_cover](frozenset(ctransforms), maxrefine))
                     # confirm cover and greedily optimize order
                     mask = numpy.ones(len(levels), dtype=bool)
                     while mask.any():

From 3097885b113d26b0baf194fa07d68e013cc0f5fd Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 24 Jun 2022 14:09:45 +0200
Subject: [PATCH 03/19] remove unused simplices property

This patch removes the Reference object's simplices property which was not used
internally.
---
 nutils/element.py | 26 --------------------------
 1 file changed, 26 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index b619ff81b..10b3d2c38 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -376,10 +376,6 @@ def getpoints(self, ischeme, degree):
             return points.SimplexBezierPoints(self.ndims, degree)
         return super().getpoints(ischeme, degree)
 
-    @property
-    def simplices(self):
-        return (transform.Identity(self.ndims), self),
-
     def get_ndofs(self, degree):
         prod = lambda start, stop: functools.reduce(operator.mul, range(start, stop), 1)
         return prod(degree+1, degree+1+self.ndims) // prod(1, self.ndims+1)
@@ -690,10 +686,6 @@ def child_refs(self):
     def inside(self, point, eps=0):
         return self.ref1.inside(point[:self.ref1.ndims], eps) and self.ref2.inside(point[self.ref1.ndims:], eps)
 
-    @property
-    def simplices(self):
-        return tuple((transform.TensorChild(trans1, trans2), TensorReference(simplex1, simplex2)) for trans1, simplex1 in self.ref1.simplices for trans2, simplex2 in self.ref2.simplices)
-
     def get_ndofs(self, degree):
         return self.ref1.get_ndofs(degree)*self.ref2.get_ndofs(degree)
 
@@ -791,12 +783,6 @@ def getpoints(self, ischeme, degree):
             return points.CoordsPoints(self.vertices[[1, 2, 4, 3, 0]])
         return points.ConePoints(self.edgeref.getpoints(ischeme, degree), self.etrans, self.tip)
 
-    @property
-    def simplices(self):
-        if self.nverts == self.ndims+1 or self.edgeref.ndims == 2 and self.edgeref.nverts == 4:  # simplices and square-based pyramids are ok
-            return [(transform.Identity(self.ndims), self)]
-        return tuple((transform.Identity(self.ndims), ref.cone(self.etrans*trans, self.tip)) for trans, ref in self.edgeref.simplices)
-
     def inside(self, point, eps=0):
         # point = etrans.apply(epoint) * xi + tip * (1-xi) => etrans.apply(epoint) = tip + (point-tip) / xi
         xi = numpy.dot(self.etrans.ext, point-self.tip) / numpy.dot(self.etrans.ext, self.etrans.offset-self.tip)
@@ -831,10 +817,6 @@ def edge_refs(self):
     def getpoints(self, ischeme, degree):
         return self.baseref.getpoints(ischeme, degree)
 
-    @property
-    def simplices(self):
-        return self.baseref.simplices
-
     def get_ndofs(self, degree):
         return self.baseref.get_ndofs(degree)
 
@@ -929,10 +911,6 @@ def getpoints(self, ischeme, degree):
         childpoints = [points.TransformPoints(ref.getpoints(ischeme, degree), trans) for trans, ref in self.children if ref]
         return points.ConcatPoints(childpoints, points.find_duplicates(childpoints) if dedup else ())
 
-    @property
-    def simplices(self):
-        return [(trans2*trans1, simplex) for trans2, child in self.children for trans1, simplex in (child.simplices if child else [])]
-
     @property
     def edge_transforms(self):
         return tuple(self.baseref.edge_transforms) \
@@ -1069,10 +1047,6 @@ def _nlinear_by_level(self, n):
     def subrefs(self):
         return tuple(ref.cone(trans, self._midpoint) for trans, ref in zip(self.baseref.edge_transforms, self._edge_refs) if ref)
 
-    @property
-    def simplices(self):
-        return [simplex for subvol in self.subrefs for simplex in subvol.simplices]
-
     def getpoints(self, ischeme, degree):
         if ischeme == 'vertex':
             return self.baseref.getpoints(ischeme, degree)

From cea14b7f220dddf064c243f81450d13f982c89d2 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 4 Jul 2022 08:15:14 +0200
Subject: [PATCH 04/19] remove unnessecary volume check

This patch removes an unnecessary consistency check in the Reference object's
slice method.
---
 nutils/element.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 10b3d2c38..2d3c1b6d0 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -226,8 +226,7 @@ def slice(self, levelfunc, ndivisions):
         if not any(refs):
             return self.empty
 
-        mosaic = MosaicReference(self, refs, midpoint)
-        return self.empty if mosaic.volume == 0 else mosaic if mosaic.volume < self.volume else self
+        return MosaicReference(self, refs, midpoint)
 
     def cone(self, trans, tip):
         return Cone(self, trans, tip)

From 064418a6795d18420f29869aed8b52b438d418f4 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 14 Nov 2022 11:52:53 +0100
Subject: [PATCH 05/19] enhance test for non-destructive slice

A slice operation returns the empty element in case only a single edge is
changed, or, conversely, the full element in case all but one edges are
unchanged. This patch relaxes this criterion for elements of 3D and up by
extending it to ndims-1 edges. The reasons for this are twofold: 1. the test is
more powerful, resulting in fewer unnecessary mosaic elements, and 2. the new
criterion generalizes to 1D, which is used in a subsequent commit to simplify
the slice method.
---
 nutils/element.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 2d3c1b6d0..60735a43d 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -207,9 +207,9 @@ def slice(self, levelfunc, ndivisions):
         else:
 
             refs = [edgeref.slice(lambda vertices: levelfunc(edgetrans.apply(vertices)), ndivisions) for edgetrans, edgeref in self.edges]
-            if sum(ref != baseref for ref, baseref in zip(refs, self.edge_refs)) <= 1:
+            if sum(ref != baseref for ref, baseref in zip(refs, self.edge_refs)) < self.ndims:
                 return self
-            if sum(bool(ref) for ref in refs) <= 1:
+            if sum(bool(ref) for ref in refs) < self.ndims:
                 return self.empty
 
             clevel = levelfunc(self.centroid[_])[0]

From 7b354b70ba8b326af53aae17b736e5952cae8bf0 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 4 Jul 2022 07:12:43 +0200
Subject: [PATCH 06/19] simplify slice logic

This patch extends the Reference object's slice method's recursion to the level
of points, which centralizes the check for corner cases in which an element is
entirely preserved or entirely removed.
---
 nutils/element.py | 21 +++++++--------------
 1 file changed, 7 insertions(+), 14 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 60735a43d..171fb6440 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -191,6 +191,13 @@ def slice(self, levelfunc, ndivisions):
             return self
         if numpy.less_equal(levels, 0).all():
             return self.empty
+        assert self.ndims >= 1
+
+        refs = [edgeref.slice(lambda vertices: levelfunc(edgetrans.apply(vertices)), ndivisions) for edgetrans, edgeref in self.edges]
+        if sum(ref != baseref for ref, baseref in zip(refs, self.edge_refs)) < self.ndims:
+            return self
+        if sum(bool(ref) for ref in refs) < self.ndims:
+            return self.empty
 
         nbins = 2**ndivisions
 
@@ -202,30 +209,16 @@ def slice(self, levelfunc, ndivisions):
                 return self.empty if xi == 0 and l1 < 0 or xi == nbins and l0 < 0 else self
             v0, v1 = self.vertices
             midpoint = v0 + (xi/nbins) * (v1-v0)
-            refs = [edgeref if levelfunc(edgetrans.apply(numpy.zeros((1, 0)))) > 0 else edgeref.empty for edgetrans, edgeref in self.edges]
 
         else:
 
-            refs = [edgeref.slice(lambda vertices: levelfunc(edgetrans.apply(vertices)), ndivisions) for edgetrans, edgeref in self.edges]
-            if sum(ref != baseref for ref, baseref in zip(refs, self.edge_refs)) < self.ndims:
-                return self
-            if sum(bool(ref) for ref in refs) < self.ndims:
-                return self.empty
-
             clevel = levelfunc(self.centroid[_])[0]
-
             select = clevel*levels <= 0 if clevel != 0 else levels != 0
             levels = levels[select]
             vertices = self.vertices[select]
-
             xi = numpy.round(levels/(levels-clevel) * nbins)
             midpoint = numpy.mean(vertices + (self.centroid-vertices)*(xi/nbins)[:, _], axis=0)
 
-        if tuple(refs) == tuple(self.edge_refs):
-            return self
-        if not any(refs):
-            return self.empty
-
         return MosaicReference(self, refs, midpoint)
 
     def cone(self, trans, tip):

From 2ddc5c1954b3f585ffdcaf75cf4442e55c21ade8 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 7 Jul 2022 16:26:15 +0200
Subject: [PATCH 07/19] improve MosaicReference.vertices and move to init

This patch makes the MosaicReference object's vertices attribute more efficient
by relying on the property that the new edges have their first vertices in
common with those of the original ones. Furthermore, attribute construction is
moved to the constructor, in anticipation of a future commit that will add
simultaneous tracking of vertex-edge relations.
---
 nutils/element.py | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 171fb6440..cde0d2b84 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -944,16 +944,19 @@ def get_edge_dofs(self, degree, iedge):
 class MosaicReference(Reference):
     'triangulation'
 
-    __slots__ = 'baseref', '_edge_refs', '_midpoint', 'edge_refs', 'edge_transforms'
-    __cache__ = 'vertices', 'subrefs'
+    __slots__ = 'baseref', '_edge_refs', '_midpoint', 'edge_refs', 'edge_transforms', 'vertices'
+    __cache__ = 'subrefs'
 
     @types.apply_annotations
     def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
         assert len(edge_refs) == baseref.nedges
         assert edge_refs != tuple(baseref.edge_refs)
+        assert midpoint.shape == (baseref.ndims,)
+
+        vertices = list(baseref.vertices)
+        assert not any((baseref.vertices == midpoint).all(1))
+        vertices.append(numpy.asarray(midpoint))
 
-        self.baseref = baseref
-        self._edge_refs = edge_refs
         self._midpoint = numpy.asarray(midpoint)
         self.edge_refs = list(edge_refs)
         self.edge_transforms = list(baseref.edge_transforms)
@@ -967,6 +970,11 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
 
         else:
 
+            for trans, edge, newedge in zip(baseref.edge_transforms, baseref.edge_refs, edge_refs):
+                for v in trans.apply(newedge.vertices[edge.nverts:]):
+                    if not any((v == v_).all() for v_ in vertices[baseref.nverts:]):
+                        vertices.append(v)
+
             newedges = [(etrans1, etrans2, edge) for (etrans1, orig), new in zip(baseref.edges, edge_refs) for etrans2, edge in new.edges[orig.nedges:]]
             for (iedge1, iedge2), (jedge1, jedge2) in baseref.ribbons:
                 Ei = edge_refs[iedge1]
@@ -989,12 +997,11 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
                 self.edge_transforms.append(newtrans)
                 self.edge_refs.append(edge.cone(extrudetrans, tip))
 
-        super().__init__(baseref.ndims)
+        self.baseref = baseref
+        self._edge_refs = edge_refs
+        self.vertices = types.frozenarray(vertices, copy=False)
 
-    @property
-    def vertices(self):
-        vertices, indices = util.unique([vertex for etrans, eref in self.edges if eref for vertex in etrans.apply(eref.vertices)], key=types.arraydata)
-        return types.frozenarray(vertices)
+        super().__init__(baseref.ndims)
 
     def __and__(self, other):
         if other in (self, self.baseref):

From 7c879cd1145fe47cfdb859fe57ae3ec1522e6180 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 7 Jul 2022 12:44:13 +0200
Subject: [PATCH 08/19] add edge_vertices

This commit add the Reference object's edge_vertices attribute, which is a
tuple of integer arrays identifying edge vertices with those of the element.
---
 nutils/element.py     | 46 ++++++++++++++++++++++++++++++++++++++++---
 tests/test_element.py |  5 +++++
 2 files changed, 48 insertions(+), 3 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index cde0d2b84..e1a1996d8 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -35,6 +35,30 @@ def __init__(self, ndims: int):
     def nverts(self):
         return len(self.vertices)
 
+    @property
+    def vertices(self):
+        raise NotImplementedError(self)
+
+    @property
+    def edge_vertices(self):
+        '''Relation between edge and volume vertices.
+
+        Given a volume reference `vol`, `vol.edge_vertices` is a sequence of
+        integer arrays that specifies per edge (outer sequence, corresponding
+        to `vol.edges`) for each vertex (inner sequence, corresponding to
+        `vol.edges[iedge].vertices`) its index in `vol` (corresponding to
+        `vol.vertices`).
+        '''
+
+        # naive implementation; will be removed in a later commit
+        edge_vertices = []
+        for etrans, eref in self.edges:
+            dist = numpy.linalg.norm(self.vertices[:,_,:] - etrans.apply(eref.vertices), axis=2)
+            emap = numpy.argmin(dist, 0)
+            assert (dist[emap, numpy.arange(eref.nverts)] < 1e-15).all(), dist
+            edge_vertices.append(types.frozenarray(emap, copy=False))
+        return tuple(edge_vertices)
+
     __and__ = lambda self, other: self if self == other else NotImplemented
     __or__ = lambda self, other: self if self == other else NotImplemented
     __rand__ = lambda self, other: self.__and__(other)
@@ -298,6 +322,10 @@ def __init__(self, baseref: strictreference):
     def vertices(self):
         return self.baseref.vertices
 
+    @property
+    def edge_vertices(self):
+        return self.baseref.edge_vertices
+
     @property
     def edge_transforms(self):
         return self.baseref.edge_transforms
@@ -329,11 +357,15 @@ class SimplexReference(Reference):
     'simplex reference'
 
     __slots__ = ()
-    __cache__ = 'edge_refs', 'edge_transforms', 'ribbons', '_get_poly_coeffs_bernstein', '_get_poly_coeffs_lagrange', '_integer_barycentric_coordinates'
+    __cache__ = 'edge_refs', 'edge_transforms', 'edge_vertices', 'ribbons', '_get_poly_coeffs_bernstein', '_get_poly_coeffs_lagrange', '_integer_barycentric_coordinates'
 
     @property
     def vertices(self):
-        return types.frozenarray(numpy.concatenate([numpy.zeros(self.ndims)[_, :], numpy.eye(self.ndims)], axis=0), copy=False)
+        return types.frozenarray(numpy.eye(self.ndims+1)[1:].T, copy=False) # first vertex in origin
+
+    @property
+    def edge_vertices(self):
+        return tuple(types.frozenarray(numpy.arange(self.ndims+1).repeat(self.ndims).reshape(self.ndims,self.ndims+1).T[::-1], copy=False))
 
     @property
     def edge_refs(self):
@@ -564,7 +596,7 @@ class TensorReference(Reference):
     'tensor reference'
 
     __slots__ = 'ref1', 'ref2'
-    __cache__ = 'vertices', 'edge_transforms', 'ribbons', 'child_transforms', 'getpoints', 'get_poly_coeffs', 'centroid'
+    __cache__ = 'vertices', 'edge_transforms', 'edge_vertices', 'ribbons', 'child_transforms', 'getpoints', 'get_poly_coeffs', 'centroid'
 
     def __init__(self, ref1, ref2):
         assert not isinstance(ref1, TensorReference)
@@ -582,6 +614,14 @@ def vertices(self):
         vertices[:, :, self.ref1.ndims:] = self.ref2.vertices[_, :]
         return types.frozenarray(vertices.reshape(self.ref1.nverts*self.ref2.nverts, self.ndims), copy=False)
 
+    @property
+    def edge_vertices(self):
+        n1 = self.ref1.nverts
+        n2 = self.ref2.nverts
+        edge_vertices = [everts[:,_] * n2 + numpy.arange(n2) for everts in self.ref1.edge_vertices] \
+                      + [numpy.arange(n1)[:,_] * n2 + everts for everts in self.ref2.edge_vertices]
+        return tuple(types.frozenarray(e.ravel(), copy=False) for e in edge_vertices)
+
     @property
     def centroid(self):
         return types.frozenarray(numpy.concatenate([self.ref1.centroid, self.ref2.centroid]), copy=False)
diff --git a/tests/test_element.py b/tests/test_element.py
index 8d5ee6735..8271908ee 100644
--- a/tests/test_element.py
+++ b/tests/test_element.py
@@ -72,6 +72,11 @@ def test_inside(self):
         if self.ref.ndims:
             self.assertFalse(self.ref.inside(-numpy.ones(self.ref.ndims)))
 
+    @parametrize.enable_if(lambda ref, **kwargs: not isinstance(ref, element.WithChildrenReference) and ref.ndims >= 1)
+    def test_edge_vertices(self):
+        for etrans, eref, everts in zip(self.ref.edge_transforms, self.ref.edge_refs, self.ref.edge_vertices):
+            self.assertAllEqual(self.ref.vertices[everts], etrans.apply(eref.vertices))
+
 
 elem('point', ref=element.PointReference(), exactcentroid=numpy.zeros((0,)))
 elem('point2', ref=element.PointReference()**2, exactcentroid=numpy.zeros((0,)))

From b7ec97238fcad8deb3631c4736fed30845455b5d Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 4 Jul 2022 07:24:21 +0200
Subject: [PATCH 09/19] fix 1D multitrim

This patch accounts for the situation that a 1D element may have more than two
vertices. This is the case when a line is trimmed for a second time: the mosaic
resulting from the first slice has three edges, one of them empty. A unittest
is added to cover this situation.
---
 nutils/element.py        |  5 +++--
 tests/test_finitecell.py | 10 +++++++++-
 2 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index e1a1996d8..fa330547f 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -227,11 +227,12 @@ def slice(self, levelfunc, ndivisions):
 
         if self.ndims == 1:
 
-            l0, l1 = levels
+            iedge = [i for (i,), edge in zip(self.edge_vertices, self.edge_refs) if edge]
+            l0, l1 = levels[iedge]
             xi = numpy.round(l0/(l0-l1) * nbins)
             if xi in (0, nbins):
                 return self.empty if xi == 0 and l1 < 0 or xi == nbins and l0 < 0 else self
-            v0, v1 = self.vertices
+            v0, v1 = self.vertices[iedge]
             midpoint = v0 + (xi/nbins) * (v1-v0)
 
         else:
diff --git a/tests/test_finitecell.py b/tests/test_finitecell.py
index 51a192e48..5e103af1a 100644
--- a/tests/test_finitecell.py
+++ b/tests/test_finitecell.py
@@ -260,7 +260,15 @@ def test_locate(self):
 
 class multitrim(TestCase):
 
-    def test(self):
+    def test_1d(self):
+        domain, geom = mesh.rectilinear([3])
+        trimmed = domain.trim(geom-1.2, maxrefine=0).trim(1.8-geom, maxrefine=0)
+        self.assertEqual(len(trimmed), 1) # trimmed consists of a single line mosaic with two new edges
+        trimmed.check_boundary(geom, elemwise=True, print=self.fail)
+        L = trimmed.integrate(function.J(geom), ischeme='gauss1')
+        numpy.testing.assert_almost_equal(L, .6, decimal=3)
+
+    def test_2d(self):
         domain, geom = mesh.rectilinear([[-1, 1], [-1, 1]])
         geom_rel = (function.rotmat(numpy.pi/6) * geom).sum(-1)
         for itrim in range(4):

From d92102f1f6fd96fa947c3611564d09bca9d15c00 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Thu, 7 Jul 2022 17:01:46 +0200
Subject: [PATCH 10/19] add simplices

This patch adds the Reference object's simplices attribute, which partitions
an element in a sequence of simplices referencing its vertices.
---
 nutils/element.py     | 60 ++++++++++++++++++++++++++++++++++++++++---
 tests/test_element.py | 14 +++++++++-
 2 files changed, 70 insertions(+), 4 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index fa330547f..ff139696e 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -37,6 +37,18 @@ def nverts(self):
 
     @property
     def vertices(self):
+        '''Vertex coordinates.'''
+
+        raise NotImplementedError(self)
+
+    @property
+    def simplices(self):
+        '''Partition of the element consisting of simplices.
+
+        The `simplices` attribute is a sequence of integer arrays that specify
+        per simplex the indices of the vertices in :attr:`vertices`.
+        '''
+
         raise NotImplementedError(self)
 
     @property
@@ -327,6 +339,10 @@ def vertices(self):
     def edge_vertices(self):
         return self.baseref.edge_vertices
 
+    @property
+    def simplices(self):
+        return types.frozenarray(numpy.empty([0, self.ndims+1], dtype=int), copy=False)
+
     @property
     def edge_transforms(self):
         return self.baseref.edge_transforms
@@ -368,6 +384,10 @@ def vertices(self):
     def edge_vertices(self):
         return tuple(types.frozenarray(numpy.arange(self.ndims+1).repeat(self.ndims).reshape(self.ndims,self.ndims+1).T[::-1], copy=False))
 
+    @property
+    def simplices(self):
+        return types.frozenarray(numpy.arange(self.ndims+1)[numpy.newaxis], copy=False)
+
     @property
     def edge_refs(self):
         assert self.ndims > 0
@@ -597,7 +617,7 @@ class TensorReference(Reference):
     'tensor reference'
 
     __slots__ = 'ref1', 'ref2'
-    __cache__ = 'vertices', 'edge_transforms', 'edge_vertices', 'ribbons', 'child_transforms', 'getpoints', 'get_poly_coeffs', 'centroid'
+    __cache__ = 'vertices', 'edge_transforms', 'edge_vertices', 'ribbons', 'child_transforms', 'getpoints', 'get_poly_coeffs', 'centroid', 'simplices'
 
     def __init__(self, ref1, ref2):
         assert not isinstance(ref1, TensorReference)
@@ -623,6 +643,25 @@ def edge_vertices(self):
                       + [numpy.arange(n1)[:,_] * n2 + everts for everts in self.ref2.edge_vertices]
         return tuple(types.frozenarray(e.ravel(), copy=False) for e in edge_vertices)
 
+    @property
+    def simplices(self):
+        if self.ref1.ndims != 1 and self.ref2.ndims != 1:
+            raise NotImplementedError((self.ref1, self.ref2))
+        # For an n-dimensional simplex with vertices a0,a1,..,an, the extruded
+        # element has vertices a0,a1,..,an,b0,b1,..,bn. These can be divided in
+        # simplices by selecting a0,a1,..,an,b0; a1,..,an,b0,n1; and so on until
+        # an,b0,b1,..,bn; resulting in n+1 n+1-dimensional simplices.
+        indices = self.ref1.simplices[:,numpy.newaxis,:,numpy.newaxis] * self.ref2.nverts \
+                + self.ref2.simplices[numpy.newaxis,:,numpy.newaxis,:]
+        if self.ref1.ndims != 1:
+            indices = indices.swapaxes(2,3) # simplex strips require penultimate axis to be of length 2
+        assert indices.shape[2] == 2
+        indices = numeric.overlapping(indices.reshape(-1, 2*self.ndims), n=self.ndims+1).copy() # nsimplex x nstrip x ndims+1
+        # to see determinants: X = self.vertices[indices]; numpy.linalg.det(X[...,1:,:] - X[...,:1,:])
+        if self.ndims % 2 == 0: # simplex strips of even dimension (e.g. triangles) have alternating orientation
+            indices[:,::2,:2] = indices[:,::2,1::-1].copy() # flip every other simplex
+        return types.frozenarray(indices.reshape(-1, self.ndims+1), copy=False)
+
     @property
     def centroid(self):
         return types.frozenarray(numpy.concatenate([self.ref1.centroid, self.ref2.centroid]), copy=False)
@@ -769,6 +808,10 @@ def __init__(self, edgeref, etrans, tip: types.arraydata):
     def vertices(self):
         return types.frozenarray(numpy.vstack([[self.tip], self.etrans.apply(self.edgeref.vertices)]), copy=False)
 
+    @property
+    def simplices(self):
+        return types.frozenarray([[0, *emap[simplex]] for emap, eref in zip(self.edge_vertices, self.edge_refs) for simplex in eref.simplices])
+
     @property
     def edge_transforms(self):
         edge_transforms = [self.etrans]
@@ -985,8 +1028,8 @@ def get_edge_dofs(self, degree, iedge):
 class MosaicReference(Reference):
     'triangulation'
 
-    __slots__ = 'baseref', '_edge_refs', '_midpoint', 'edge_refs', 'edge_transforms', 'vertices'
-    __cache__ = 'subrefs'
+    __slots__ = 'baseref', '_edge_refs', '_midpoint', 'edge_refs', 'edge_transforms', 'vertices', '_imidpoint'
+    __cache__ = 'subrefs', 'simplices'
 
     @types.apply_annotations
     def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
@@ -997,6 +1040,7 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
         vertices = list(baseref.vertices)
         assert not any((baseref.vertices == midpoint).all(1))
         vertices.append(numpy.asarray(midpoint))
+        imidpoint = baseref.nverts
 
         self._midpoint = numpy.asarray(midpoint)
         self.edge_refs = list(edge_refs)
@@ -1041,6 +1085,7 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
         self.baseref = baseref
         self._edge_refs = edge_refs
         self.vertices = types.frozenarray(vertices, copy=False)
+        self._imidpoint = imidpoint
 
         super().__init__(baseref.ndims)
 
@@ -1083,6 +1128,15 @@ def __rsub__(self, other):
     def _nlinear_by_level(self, n):
         return self.baseref._nlinear_by_level(n)
 
+    @property
+    def simplices(self):
+        indices = []
+        for vmap, etrans, eref in zip(self.edge_vertices, self.baseref.edge_transforms, self._edge_refs):
+            for index in vmap[eref.simplices]:
+                assert self._imidpoint not in index
+                indices.append([self._imidpoint, *index] if not etrans.isflipped else [index[0], self._imidpoint, *index[1:]])
+        return types.frozenarray(indices, dtype=int)
+
     @property
     def subrefs(self):
         return tuple(ref.cone(trans, self._midpoint) for trans, ref in zip(self.baseref.edge_transforms, self._edge_refs) if ref)
diff --git a/tests/test_element.py b/tests/test_element.py
index 8271908ee..dabca0098 100644
--- a/tests/test_element.py
+++ b/tests/test_element.py
@@ -1,6 +1,6 @@
 from nutils import element, _util as util
 from nutils.testing import TestCase, parametrize
-import numpy
+import numpy, math
 
 
 @parametrize
@@ -77,6 +77,18 @@ def test_edge_vertices(self):
         for etrans, eref, everts in zip(self.ref.edge_transforms, self.ref.edge_refs, self.ref.edge_vertices):
             self.assertAllEqual(self.ref.vertices[everts], etrans.apply(eref.vertices))
 
+    @parametrize.enable_if(lambda ref, **kwargs: not isinstance(ref, element.WithChildrenReference) and ref.ndims >= 1)
+    def test_simplices(self):
+        volume = 0
+        centroid = 0
+        for simplex in self.ref.vertices[self.ref.simplices]:
+            simplex_volume = numpy.linalg.det(simplex[1:] - simplex[0]) / math.factorial(self.ref.ndims)
+            self.assertGreater(simplex_volume, 0)
+            volume += simplex_volume
+            centroid += simplex.mean(axis=0) * simplex_volume
+        centroid /= volume
+        self.assertAlmostEqual(volume, self.ref.volume)
+        self.assertAllAlmostEqual(centroid, self.exactcentroid)
 
 elem('point', ref=element.PointReference(), exactcentroid=numpy.zeros((0,)))
 elem('point2', ref=element.PointReference()**2, exactcentroid=numpy.zeros((0,)))

From cfa9d09af9968a4a0cb33639592406574f3bbf87 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 8 Jul 2022 09:57:51 +0200
Subject: [PATCH 11/19] add simplex_transforms

This patch adds the Reference object's simplex_transforms attribute, which
converts the simplices into transform items to be combined with a simplex
reference.
---
 nutils/element.py     | 19 +++++++++++++++++++
 tests/test_element.py |  6 ++++++
 2 files changed, 25 insertions(+)

diff --git a/nutils/element.py b/nutils/element.py
index ff139696e..6ef148dff 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -51,6 +51,19 @@ def simplices(self):
 
         raise NotImplementedError(self)
 
+    @property
+    def simplex_transforms(self):
+        '''Sequence of transforms from simplex to parent element.
+
+        The `simplex_transforms` attribute is a sequence of objects of type
+        :class:`nutils.transform.TransformItem` that provide per simplex the
+        coordinate mapping from the simplex to the parent element. The origin
+        of the simplex-local coordinate system maps to its first vertex, the
+        first unit vector to the second, the second to the third, and so on.
+        '''
+
+        return tuple(transform.Square((vertices[1:] - vertices[0]).T, vertices[0]) for vertices in self.vertices[self.simplices])
+
     @property
     def edge_vertices(self):
         '''Relation between edge and volume vertices.
@@ -388,6 +401,12 @@ def edge_vertices(self):
     def simplices(self):
         return types.frozenarray(numpy.arange(self.ndims+1)[numpy.newaxis], copy=False)
 
+    @property
+    def simplex_transforms(self):
+        # The definition of self.vertices is such that the conventions of
+        # Reference.simplex_transforms result in the identity map.
+        return transform.Identity(self.ndims),
+
     @property
     def edge_refs(self):
         assert self.ndims > 0
diff --git a/tests/test_element.py b/tests/test_element.py
index dabca0098..f389dd250 100644
--- a/tests/test_element.py
+++ b/tests/test_element.py
@@ -90,6 +90,12 @@ def test_simplices(self):
         self.assertAlmostEqual(volume, self.ref.volume)
         self.assertAllAlmostEqual(centroid, self.exactcentroid)
 
+    @parametrize.enable_if(lambda ref, **kwargs: not isinstance(ref, element.WithChildrenReference) and ref.ndims >= 1)
+    def test_simplex_transforms(self):
+        for simplex, strans in zip(self.ref.vertices[self.ref.simplices], self.ref.simplex_transforms):
+            self.assertAllEqual(strans.linear, (simplex[1:] - simplex[0]).T)
+            self.assertAllEqual(strans.offset, simplex[0])
+
 elem('point', ref=element.PointReference(), exactcentroid=numpy.zeros((0,)))
 elem('point2', ref=element.PointReference()**2, exactcentroid=numpy.zeros((0,)))
 elem('line', ref=element.LineReference(), exactcentroid=numpy.array([.5]))

From d957b3a0147f0e15e4dc24118dc3a2d238ca68c0 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 8 Jul 2022 10:27:19 +0200
Subject: [PATCH 12/19] implement inside based in simplex_transforms

This patch adds a base implementation of Reference.inside using simplices,
replacing the implementations of Empty, Simplex, and Mosaic.
---
 nutils/element.py | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 6ef148dff..419446d5f 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -64,6 +64,14 @@ def simplex_transforms(self):
 
         return tuple(transform.Square((vertices[1:] - vertices[0]).T, vertices[0]) for vertices in self.vertices[self.simplices])
 
+    def inside(self, point, eps=0):
+        for strans in self.simplex_transforms:
+            spoint = strans.invapply(point) # point in simplex coordinates
+            tol = -eps / strans.det # account for simplex scale
+            if all(bary >= tol for bary in (*spoint, 1-spoint.sum())):
+                return True
+        return False
+
     @property
     def edge_vertices(self):
         '''Relation between edge and volume vertices.
@@ -379,9 +387,6 @@ def child_refs(self):
     def trim(self, levels, maxrefine, ndivisions):
         return self
 
-    def inside(self, point, eps=0):
-        return False
-
 
 class SimplexReference(Reference):
     'simplex reference'
@@ -484,9 +489,6 @@ def _get_poly_coeffs_lagrange(self, degree):
     def get_edge_dofs(self, degree, iedge):
         return types.frozenarray(tuple(i for i, j in enumerate(self._integer_barycentric_coordinates(degree)) if j[iedge] == 0), dtype=int)
 
-    def inside(self, point, eps=0):
-        return numpy.greater_equal(point, -eps).all(axis=0) and numpy.less_equal(numpy.sum(point, axis=0), 1+eps)
-
 
 class PointReference(SimplexReference):
     '0D simplex'
@@ -1173,9 +1175,6 @@ def getpoints(self, ischeme, degree):
         # later, to be combined with a reduction of gauss schemes of any degree.
         return points.ConcatPoints(subpoints, dups)
 
-    def inside(self, point, eps=0):
-        return any(subref.inside(point, eps=eps) for subref in self.subrefs)
-
     def get_ndofs(self, degree):
         return self.baseref.get_ndofs(degree)
 

From 198697913685e483c36b926b604d22af8eb9cdd7 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 8 Jul 2022 11:31:37 +0200
Subject: [PATCH 13/19] remove subrefs

This patch removes the MosaicReference object's subrefs attribute and rewrites
its points method in terms of simplices.
---
 nutils/element.py | 20 +++++++-------------
 1 file changed, 7 insertions(+), 13 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 419446d5f..6079534fb 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -1050,7 +1050,7 @@ class MosaicReference(Reference):
     'triangulation'
 
     __slots__ = 'baseref', '_edge_refs', '_midpoint', 'edge_refs', 'edge_transforms', 'vertices', '_imidpoint'
-    __cache__ = 'subrefs', 'simplices'
+    __cache__ = 'simplices'
 
     @types.apply_annotations
     def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
@@ -1158,22 +1158,16 @@ def simplices(self):
                 indices.append([self._imidpoint, *index] if not etrans.isflipped else [index[0], self._imidpoint, *index[1:]])
         return types.frozenarray(indices, dtype=int)
 
-    @property
-    def subrefs(self):
-        return tuple(ref.cone(trans, self._midpoint) for trans, ref in zip(self.baseref.edge_transforms, self._edge_refs) if ref)
-
     def getpoints(self, ischeme, degree):
         if ischeme == 'vertex':
             return self.baseref.getpoints(ischeme, degree)
-        if ischeme == '_centroid':
+        elif ischeme in ('gauss', 'uniform', 'bezier'):
+            simplexpoints = getsimplex(self.ndims).getpoints(ischeme, degree)
+            subpoints = [points.TransformPoints(simplexpoints, strans) for strans in self.simplex_transforms]
+            dups = points.find_duplicates(subpoints) if ischeme == 'bezier' else ()
+            return points.ConcatPoints(subpoints, dups)
+        else:
             return super().getpoints(ischeme, degree)
-        subpoints = [subvol.getpoints(ischeme, degree) for subvol in self.subrefs]
-        dups = points.find_duplicates(subpoints) if ischeme == 'bezier' else ()
-        # NOTE We could consider postprocessing gauss1 to a single point scheme,
-        # rather than a concatenation, but for that we would have to verify that
-        # the centroid is contained in the element. We leave this optimization for
-        # later, to be combined with a reduction of gauss schemes of any degree.
-        return points.ConcatPoints(subpoints, dups)
 
     def get_ndofs(self, degree):
         return self.baseref.get_ndofs(degree)

From 495e6dbb18d7fd12bab61d89906f097a50211788 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 8 Jul 2022 15:01:56 +0200
Subject: [PATCH 14/19] change mosaic edges from cone to simplex

This patch changes the MosaicReference object's constructor to generate trimmed
edges from simplices rather than cone objects.
---
 nutils/element.py | 143 +++++++++++++++++++++++++++++++++++-----------
 1 file changed, 111 insertions(+), 32 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 6079534fb..8d6636ab6 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -1049,7 +1049,7 @@ def get_edge_dofs(self, degree, iedge):
 class MosaicReference(Reference):
     'triangulation'
 
-    __slots__ = 'baseref', '_edge_refs', '_midpoint', 'edge_refs', 'edge_transforms', 'vertices', '_imidpoint'
+    __slots__ = 'baseref', '_edge_refs', '_midpoint', 'edge_refs', 'edge_transforms', 'vertices', '_imidpoint', 'edge_vertices'
     __cache__ = 'simplices'
 
     @types.apply_annotations
@@ -1057,56 +1057,135 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
         assert len(edge_refs) == baseref.nedges
         assert edge_refs != tuple(baseref.edge_refs)
         assert midpoint.shape == (baseref.ndims,)
+        assert not any((baseref.vertices == midpoint).all(1))
+        assert all(numpy.all(edge.vertices == newedge.vertices[:edge.nverts])
+            for edge, newedge in zip(baseref.edge_refs, edge_refs))
 
         vertices = list(baseref.vertices)
-        assert not any((baseref.vertices == midpoint).all(1))
         vertices.append(numpy.asarray(midpoint))
         imidpoint = baseref.nverts
 
-        self._midpoint = numpy.asarray(midpoint)
-        self.edge_refs = list(edge_refs)
-        self.edge_transforms = list(baseref.edge_transforms)
+        # The remainder of this constructor is concerned with establishing the
+        # edges of the mosaic, and setting up the corresponding edge-vertex
+        # relations.
 
         if baseref.ndims == 1:
 
-            assert any(edge_refs) and not all(edge_refs), 'invalid 1D mosaic: exactly one edge should be non-empty'
-            iedge, = [i for i, edge in enumerate(edge_refs) if edge]
-            self.edge_refs.append(getsimplex(0))
-            self.edge_transforms.append(transform.Updim(linear=numpy.zeros((1, 0)), offset=self._midpoint, isflipped=not baseref.edge_transforms[iedge].isflipped))
+            # For 1D elements the situation is simple: the midpoint represents
+            # the only new vertex (already added to vertices) as well as the
+            # only new edge, with a trivial edge-vertex relationship.
+
+            edge_vertices = (*baseref.edge_vertices, types.frozenarray([imidpoint]))
+            orientation = [not trans.isflipped for trans, edge in zip(baseref.edge_transforms, edge_refs) if edge]
+            assert len(orientation) == 1, 'invalid 1D mosaic: exactly one edge should be non-empty'
 
         else:
 
-            for trans, edge, newedge in zip(baseref.edge_transforms, baseref.edge_refs, edge_refs):
+            # For higher-dimensional elements the situation is more complex.
+            # Firstly, the new edge_refs may introduce new vertices. Luckily
+            # here the previously asserted convention applies that the vertices
+            # of the original edge are repeated in the modified edge, so we can
+            # focus our attention on the new ones, having only to deduplicate
+            # between them.
+
+            edge_vertices = []
+            for trans, edge, emap, newedge in zip(baseref.edge_transforms, baseref.edge_refs, baseref.edge_vertices, edge_refs):
                 for v in trans.apply(newedge.vertices[edge.nverts:]):
-                    if not any((v == v_).all() for v_ in vertices[baseref.nverts:]):
+                    for i, v_ in enumerate(vertices[baseref.nverts:], start=baseref.nverts):
+                        if (v == v_).all():
+                            break # the new vertex was already added by a previous edge
+                    else:
+                        i = len(vertices)
                         vertices.append(v)
+                    emap = types.frozenarray([*emap, i])
+                edge_vertices.append(emap)
+
+            # Secondly, new edges (which will be pulled to the midpoint) can
+            # originate either from new edge-edges, or from existing ones that
+            # find themselves without a counterpart. The former situation is
+            # trivial, following the convention that existing edge transforms
+            # are copied over in the modified edge.
+
+            assert all(edge.edge_transforms == newedge.edge_transforms[:edge.nverts]
+                for edge, newedge in zip(baseref.edge_refs, edge_refs))
+
+            # The latter, however, is more tricky. This is the situation that
+            # occurs, for instance, when two out of four edges of a square are
+            # cleanly removed, making two existing edge-edges the new exterior.
+            # Identifying these situations requires an examination of all the
+            # modified edges, that is, edge-edges in locations that pre-existed
+            # in baseref. Knowing that the edges of baseref form a watertight
+            # hull, we employ the strategy of first identifying all edge-edge
+            # counterparts, and then comparing the new references in the
+            # identified locations to see if one of the two disappeared: in
+            # this case the other reference is added to the exterior set.
+
+            # NOTE: establishing edge-edge relations could potentially be
+            # cached for reuse at the level of baseref. However, since this is
+            # the only place that the information is used and all edge pairs
+            # need to anyhow be examined for gaps, it is not clear that the
+            # gains merit the additional complexity.
+
+            orientation = []
+            seen = {}
+            for edge1, newemap1, etrans1, newedge1 in zip(baseref.edge_refs, edge_vertices, baseref.edge_transforms, edge_refs):
+                newedge1_edge = zip(newedge1.edge_vertices, newedge1.edge_transforms, newedge1.edge_refs)
+                trimmed = [] # trimmed will be populated with a subset of newedge1_edge
+                for edge2, (newemap2, etrans2, newedge2) in zip(edge1.edge_refs, newedge1_edge):
+                    if edge2: # existing non-empty edge
+
+                        # To identify matching edge-edges we map their vertices
+                        # to the numbering of the baseref for comparison. Since
+                        # matching edge-edges have must have equal references,
+                        # and by construction have matching orientation, the
+                        # vertex ordering will be consistent between them.
+
+                        key = tuple(newemap1[newemap2])
+
+                        # NOTE: there have been anecdotal reports that suggest
+                        # the assumption of matching edges may be violated, but
+                        # it is not clear in what scenario this can occur. If
+                        # the 'not seen' assertion below fails, please provide
+                        # the developers with a reproducable issue for study.
+
+                        try:
+                            newedge2_ = seen.pop(key)
+                        except KeyError:
+                            seen[key] = newedge2
+                        else: # a counterpart is found, placing newedge2 against newedge2_
+                            if not newedge2:
+                                trimmed.append((newemap2, etrans2.flipped, newedge2_))
+                            elif not newedge2_:
+                                trimmed.append((newemap2, etrans2, newedge2))
+                            elif newedge2 != newedge2_:
+                                raise NotImplementedError
+
+                # Since newedge1_edge was zipped against the shorter list of
+                # original edge1.edge_refs, what remains are the new edge-edges
+                # that can be added without further examination.
+
+                trimmed.extend(newedge1_edge)
+
+                # What remains is only to extend the edge-vertex relations and
+                # to track if the new edges are left- or right-handed.
+
+                for newemap2, etrans2, newedge2 in trimmed:
+                    for simplex in newemap1[newemap2[newedge2.simplices]]:
+                        assert imidpoint not in simplex
+                        edge_vertices.append(types.frozenarray([imidpoint, *simplex]))
+                        orientation.append(not etrans1.isflipped^etrans2.isflipped)
 
-            newedges = [(etrans1, etrans2, edge) for (etrans1, orig), new in zip(baseref.edges, edge_refs) for etrans2, edge in new.edges[orig.nedges:]]
-            for (iedge1, iedge2), (jedge1, jedge2) in baseref.ribbons:
-                Ei = edge_refs[iedge1]
-                ei = Ei.edge_refs[iedge2]
-                Ej = edge_refs[jedge1]
-                ej = Ej.edge_refs[jedge2]
-                ejsubi = ej - ei
-                if ejsubi:
-                    newedges.append((self.edge_transforms[jedge1], Ej.edge_transforms[jedge2], ejsubi))
-                eisubj = ei - ej
-                if eisubj:
-                    newedges.append((self.edge_transforms[iedge1], Ei.edge_transforms[iedge2], eisubj))
-
-            extrudetrans = transform.Updim(numpy.eye(baseref.ndims-1)[:, :-1], numpy.zeros(baseref.ndims-1), isflipped=baseref.ndims % 2 == 0)
-            tip = numpy.array([0]*(baseref.ndims-2)+[1], dtype=float)
-            for etrans, trans, edge in newedges:
-                b = etrans.apply(trans.offset)
-                A = numpy.hstack([numpy.dot(etrans.linear, trans.linear), (self._midpoint-b)[:, _]])
-                newtrans = transform.Updim(A, b, isflipped=etrans.isflipped ^ trans.isflipped ^ (baseref.ndims % 2 == 1))  # isflipped logic tested up to 3D
-                self.edge_transforms.append(newtrans)
-                self.edge_refs.append(edge.cone(extrudetrans, tip))
+            assert not seen, f'leftover unmatched edges ({seen}) indicate the edges of baseref ({baseref}) are not watertight!'
 
         self.baseref = baseref
         self._edge_refs = edge_refs
         self.vertices = types.frozenarray(vertices, copy=False)
         self._imidpoint = imidpoint
+        self._midpoint = midpoint
+        self.edge_vertices = tuple(edge_vertices)
+        self.edge_refs = edge_refs + (getsimplex(baseref.ndims-1),) * len(orientation)
+        self.edge_transforms = baseref.edge_transforms + tuple(transform.Updim((vertices[1:] - vertices[0]).T, vertices[0], isflipped)
+          for vertices, isflipped in zip(self.vertices[numpy.array(edge_vertices[baseref.nedges:])], orientation))
 
         super().__init__(baseref.ndims)
 

From a26d13fd594e11ed505abb1d27dbba54e60c437b Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 8 Jul 2022 15:23:20 +0200
Subject: [PATCH 15/19] simplify setops using _with_edges helper method

This patch deduplicates some the MosaicReference object's set operations code
via the introduction of a new _with_edges helper method.
---
 nutils/element.py | 26 ++++++++++++--------------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 8d6636ab6..3bb5a49bf 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -1189,40 +1189,38 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
 
         super().__init__(baseref.ndims)
 
+    def _with_edges(self, edge_refs):
+        edge_refs = tuple(edge_refs)
+        return self.baseref if edge_refs == self.baseref.edge_refs \
+          else self.empty if not any(edge_refs) \
+          else MosaicReference(self.baseref, edge_refs, self._midpoint)
+
     def __and__(self, other):
         if other in (self, self.baseref):
             return self
         if isinstance(other, MosaicReference) and other.baseref == self:
             return other
-        if isinstance(other, MosaicReference) and self.baseref == other.baseref and numpy.equal(other._midpoint, self._midpoint).all():
-            isect_edge_refs = [selfedge & otheredge for selfedge, otheredge in zip(self._edge_refs, other._edge_refs)]
-            if not any(isect_edge_refs):
-                return self.empty
-            return MosaicReference(self.baseref, isect_edge_refs, self._midpoint)
+        if isinstance(other, MosaicReference) and self.baseref == other.baseref and other._midpoint == self._midpoint:
+            return self._with_edges(selfedge & otheredge for selfedge, otheredge in zip(self._edge_refs, other._edge_refs))
         return NotImplemented
 
     def __or__(self, other):
         if other in (self, self.baseref):
             return other
-        if isinstance(other, MosaicReference) and self.baseref == other.baseref and numpy.equal(other._midpoint, self._midpoint).all():
-            union_edge_refs = [selfedge | otheredge for selfedge, otheredge in zip(self._edge_refs, other._edge_refs)]
-            if tuple(union_edge_refs) == tuple(self.baseref.edge_refs):
-                return self.baseref
-            return MosaicReference(self.baseref, union_edge_refs, self._midpoint)
+        if isinstance(other, MosaicReference) and self.baseref == other.baseref and other._midpoint == self._midpoint:
+            return self._with_edges(selfedge | otheredge for selfedge, otheredge in zip(self._edge_refs, other._edge_refs))
         return NotImplemented
 
     def __sub__(self, other):
         if other in (self, self.baseref):
             return self.empty
         if isinstance(other, MosaicReference) and other.baseref == self:
-            inv_edge_refs = [baseedge - edge for baseedge, edge in zip(self.edge_refs, other._edge_refs)]
-            return MosaicReference(self, inv_edge_refs, other._midpoint)
+            return other._with_edges(baseedge - edge for baseedge, edge in zip(self.edge_refs, other._edge_refs))
         return NotImplemented
 
     def __rsub__(self, other):
         if other == self.baseref:
-            inv_edge_refs = [baseedge - edge for baseedge, edge in zip(other.edge_refs, self._edge_refs)]
-            return MosaicReference(other, inv_edge_refs, self._midpoint)
+            return self._with_edges(baseedge - edge for baseedge, edge in zip(other.edge_refs, self._edge_refs))
         return NotImplemented
 
     def _nlinear_by_level(self, n):

From b07e070aeae6157a2ce5141c0241552124a2b433 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 8 Jul 2022 15:31:36 +0200
Subject: [PATCH 16/19] support vertex-coinciding midpoint

This patch generalizes the MosaicReference object's implementation to correctly
support the situation that the provided midpoint coincides with a vertex of
baseref.
---
 nutils/element.py | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 3bb5a49bf..8491fa60b 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -1057,13 +1057,16 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
         assert len(edge_refs) == baseref.nedges
         assert edge_refs != tuple(baseref.edge_refs)
         assert midpoint.shape == (baseref.ndims,)
-        assert not any((baseref.vertices == midpoint).all(1))
         assert all(numpy.all(edge.vertices == newedge.vertices[:edge.nverts])
             for edge, newedge in zip(baseref.edge_refs, edge_refs))
 
         vertices = list(baseref.vertices)
-        vertices.append(numpy.asarray(midpoint))
-        imidpoint = baseref.nverts
+        match, = (baseref.vertices == midpoint).all(1).nonzero()
+        if match.size:
+            imidpoint, = match
+        else:
+            imidpoint = len(vertices)
+            vertices.append(numpy.asarray(midpoint))
 
         # The remainder of this constructor is concerned with establishing the
         # edges of the mosaic, and setting up the corresponding edge-vertex
@@ -1075,6 +1078,7 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
             # the only new vertex (already added to vertices) as well as the
             # only new edge, with a trivial edge-vertex relationship.
 
+            assert not match.size, '1D mosaic must introduce a new vertex'
             edge_vertices = (*baseref.edge_vertices, types.frozenarray([imidpoint]))
             orientation = [not trans.isflipped for trans, edge in zip(baseref.edge_transforms, edge_refs) if edge]
             assert len(orientation) == 1, 'invalid 1D mosaic: exactly one edge should be non-empty'
@@ -1171,9 +1175,9 @@ def __init__(self, baseref, edge_refs: tuple, midpoint: types.arraydata):
 
                 for newemap2, etrans2, newedge2 in trimmed:
                     for simplex in newemap1[newemap2[newedge2.simplices]]:
-                        assert imidpoint not in simplex
-                        edge_vertices.append(types.frozenarray([imidpoint, *simplex]))
-                        orientation.append(not etrans1.isflipped^etrans2.isflipped)
+                        if imidpoint not in simplex:
+                            edge_vertices.append(types.frozenarray([imidpoint, *simplex]))
+                            orientation.append(not etrans1.isflipped^etrans2.isflipped)
 
             assert not seen, f'leftover unmatched edges ({seen}) indicate the edges of baseref ({baseref}) are not watertight!'
 
@@ -1231,8 +1235,8 @@ def simplices(self):
         indices = []
         for vmap, etrans, eref in zip(self.edge_vertices, self.baseref.edge_transforms, self._edge_refs):
             for index in vmap[eref.simplices]:
-                assert self._imidpoint not in index
-                indices.append([self._imidpoint, *index] if not etrans.isflipped else [index[0], self._imidpoint, *index[1:]])
+                if self._imidpoint not in index:
+                    indices.append([self._imidpoint, *index] if not etrans.isflipped else [index[0], self._imidpoint, *index[1:]])
         return types.frozenarray(indices, dtype=int)
 
     def getpoints(self, ischeme, degree):

From e1887b3e807fa96b6577ec2da5a9951f2b973288 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 8 Jul 2022 15:28:23 +0200
Subject: [PATCH 17/19] simplify slice and remove midpoint for ndims>1

This patch simplifies the Reference object's slice method by changing the first
argument from a callable to a level array, making use of the edge_vertices
attribute to propagate this to the edges, and changes the midpoint logic for
elements of dimension >= 2 by selecting an existing vertex rather than
introducing a new midpoint.
---
 examples/platewithhole.py | 22 ++++++++++-----------
 nutils/element.py         | 41 +++++++++++++++++++++++++--------------
 tests/test_finitecell.py  |  8 ++++----
 tests/test_points.py      |  7 +++----
 4 files changed, 44 insertions(+), 34 deletions(-)

diff --git a/examples/platewithhole.py b/examples/platewithhole.py
index ca32daac5..0d46af59c 100644
--- a/examples/platewithhole.py
+++ b/examples/platewithhole.py
@@ -99,32 +99,32 @@ def test_spline(self):
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00033, places=5)
         with self.subTest('h1-error'):
-            self.assertAlmostEqual(err[1], .00672, places=5)
+            self.assertAlmostEqual(err[1], .00674, places=5)
         with self.subTest('constraints'):
             self.assertAlmostEqual64(cons, '''
                 eNpjaGBoYGBAxvrnGBow4X89g3NQFSjQwLAGq7i10Wus4k+NfM8fNWZgOGL89upc47WX0ozvXjAzPn1e
                 1TjnPACrACoJ''')
         with self.subTest('left-hand side'):
             self.assertAlmostEqual64(u, '''
-                eNpbZHbajIHhxzkGBhMgtgdi/XPypyRPvjFxO/PccPq5Vn2vcxr6luf+6xmcm2LMwLDQePf5c0bTzx8x
-                5D7vaTjnnIFhzbmlQPH5xhV39Y3vXlxtJHoh2EjvvLXR63MbgOIbjRdfrTXeecnUeO+Fn0Yrzj818j1/
-                FCh+xPjt1bnGay+lGd+9YGZ8+ryqcc55AK+AP/0=''')
+                eNpbb3bMjIHhxzkGBhMgtgdi/XMqp8RPvjLxOPPCcNq5Fn3Pcxr6luf+6xmcm2LMwLDQePf5c0bTzx8x
+                5DnvaTjnnIFhzbmlQPH5xgvu6hvfvbjaSPRCsJHeeWuj1+c2AMU3Gi++Wmu885Kp8d4LP41WnH9q5Hv+
+                KFD8iPHbq3ON115KM757wcz49HlV45zzAL8gQC8=''')
 
     def test_mixed(self):
         err, cons, u = main(nelems=4, etype='mixed', btype='std', degree=2, traction=.1, maxrefine=2, radius=.5, poisson=.3)
         with self.subTest('l2-error'):
             self.assertAlmostEqual(err[0], .00024, places=5)
         with self.subTest('h1-error'):
-            self.assertAlmostEqual(err[1], .00739, places=5)
+            self.assertAlmostEqual(err[1], .00740, places=5)
         with self.subTest('constraints'):
             self.assertAlmostEqual64(cons, '''
                 eNpjaGDADhlwiOEU1z8HZusbgukkg5BzRJqKFRoa1oD1HzfceA5NH9FmgKC10SuwOdONpM7DxDYa77gM
                 MueoMQPDEePzV2Hic42XXmoynnQRxvc3dryQbnz3Aoj91Mj3vJnx6fOqxjnnAQzkV94=''')
         with self.subTest('left-hand side'):
             self.assertAlmostEqual64(u, '''
-                eNoNzE8og3EcBvC3uUo5rNUOnBSK9/19n0Ic0Eo5oJBmRxcaB04kUnPgoETmT2w7LVrtMBy4auMw+35/
-                7/vaykFSFEopKTnIe/jU01PPU6FNWcQIn+Or5CBfSqCGD1uDYhi7/KbW+dma5aK65gX6Y8Po8HSzZQ7y
-                vBniHyvFV9aq17V7TK42O9kwFS9YUzxhjXIcZxLCnIzjTsfxah/BMFJotjUlZYz6xYeoPqEPKaigbKhb
-                9lOj9NGa9KgtVmqJH9UT36gcp71dEr6HaVS5GS8f46AcQ9itx739SQXdBL8dRqeTo1odox35poh2yJVh
-                apEueucsRWWPgpJFoLKPNzeHC/fU+yl48pDyMi6dCFbsBNJODNu2iawOoE4PoVdP4kH/UkZeaEDaUJQG
-                zMg/DouRUg==''')
+                eNoNzEEoQ3EcB/AXVymHtdqBkyLx3v/3LTQHtHJQKKHZ0YXMQS6sSM2BcrKMqbHTotUOw4GrthzWfr//
+                e6+nHJYUyUopKSnlHT717Vvfr0cpSWCWr/FVs1GuZdHKmb6QGMYRN9Qev1irXFUVTtAfG8agb5gtc5LX
+                zQj/WDm+s3b8bsBncosZZsNUvGEt8YI1w2lcSQRrMg9Pp/FmZ2EYOfTYmnIyR+PShLi+oA8pq5DsqxoH
+                qEvGaFdG1AErtclP6pnvVYnz/u4MVj2OZrfg53OceElE3Q482p9U0d0I2FGEnRK16SQdyjfFtEOuTFOv
+                DFGDi7QsxxSSIoIPGby7Jdy4l/5PxVeGeFu4dWLYtk+Rd5JI2SaKOoh2PYVRvYi6/qWCvNKE9KMqnViR
+                fyhZkYI=''')
diff --git a/nutils/element.py b/nutils/element.py
index 8491fa60b..61538167e 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -231,37 +231,38 @@ def trim(self, levels, maxrefine, ndivisions):
             else self.empty if numpy.less_equal(levels, 0).all() \
             else self.with_children(cref.trim(clevels, maxrefine-1, ndivisions)
                                     for cref, clevels in zip(self.child_refs, self.child_divide(levels, maxrefine))) if maxrefine > 0 \
-            else self.slice(lambda vertices: numeric.dot(numeric.poly_eval(self._linear_bernstein, vertices), levels), ndivisions)
+            else self.slice(numeric.poly_eval(self._linear_bernstein, self.vertices) @ levels, ndivisions)
 
     @property
     def _linear_bernstein(self):
         return self.get_poly_coeffs('bernstein', degree=1)
 
-    def slice(self, levelfunc, ndivisions):
+    def slice(self, levels, ndivisions):
         # slice along levelset by recursing over dimensions
 
-        levels = levelfunc(self.vertices)
-
-        assert numeric.isint(ndivisions)
-        assert len(levels) == self.nverts
+        assert len(levels) == len(self.vertices)
         if numpy.greater_equal(levels, 0).all():
             return self
         if numpy.less_equal(levels, 0).all():
             return self.empty
         assert self.ndims >= 1
 
-        refs = [edgeref.slice(lambda vertices: levelfunc(edgetrans.apply(vertices)), ndivisions) for edgetrans, edgeref in self.edges]
+        refs = [edgeref.slice(levels[edgeverts], ndivisions) for edgeverts, edgeref in zip(self.edge_vertices, self.edge_refs)]
         if sum(ref != baseref for ref, baseref in zip(refs, self.edge_refs)) < self.ndims:
             return self
         if sum(bool(ref) for ref in refs) < self.ndims:
             return self.empty
 
-        nbins = 2**ndivisions
-
         if self.ndims == 1:
 
+            # For 1D elements a midpoint is introduced through linear
+            # interpolation of the vertex levels, followed by a binning step to
+            # remove near-vertex cuts and improve robustness for topology-wide
+            # connectivity.
+
             iedge = [i for (i,), edge in zip(self.edge_vertices, self.edge_refs) if edge]
             l0, l1 = levels[iedge]
+            nbins = 2**ndivisions
             xi = numpy.round(l0/(l0-l1) * nbins)
             if xi in (0, nbins):
                 return self.empty if xi == 0 and l1 < 0 or xi == nbins and l0 < 0 else self
@@ -270,12 +271,22 @@ def slice(self, levelfunc, ndivisions):
 
         else:
 
-            clevel = levelfunc(self.centroid[_])[0]
-            select = clevel*levels <= 0 if clevel != 0 else levels != 0
-            levels = levels[select]
-            vertices = self.vertices[select]
-            xi = numpy.round(levels/(levels-clevel) * nbins)
-            midpoint = numpy.mean(vertices + (self.centroid-vertices)*(xi/nbins)[:, _], axis=0)
+            # For higher-dimensional elements, the first vertex that is newly
+            # introduced by an edge slice is selected to serve as 'midpoint'.
+            # In case no new vertices are introduced (all edges are either
+            # fully retained or fully removed) then the first vertex is
+            # selected that occurs in only one of the edges. Either situation
+            # guarantees that the selected vertex lies on the exterior hull.
+
+            for trans, edge, emap, newedge in zip(self.edge_transforms, self.edge_refs, self.edge_vertices, refs):
+                if newedge.nverts > edge.nverts:
+                    midpoint = trans.apply(newedge.vertices[edge.nverts])
+                    break
+            else:
+                count = numpy.zeros(self.nverts, dtype=int)
+                for emap, eref in zip(self.edge_vertices, refs):
+                    count[emap[eref.simplices]] += 1
+                midpoint = self.vertices[count==1][0]
 
         return MosaicReference(self, refs, midpoint)
 
diff --git a/tests/test_finitecell.py b/tests/test_finitecell.py
index 5e103af1a..af484829c 100644
--- a/tests/test_finitecell.py
+++ b/tests/test_finitecell.py
@@ -388,10 +388,10 @@ def test_topos(self):
         self.assertEqual(len(self.topoB), 2)
 
     def test_boundaries(self):
-        self.assertEqual(len(self.topoA.boundary), 11)
-        self.assertEqual(len(self.topoB.boundary), 8)
-        self.assertEqual(len(self.topoA.boundary['trimmed']), 5)
-        self.assertEqual(len(self.topoB.boundary['trimmed']), 5)
+        self.assertEqual(len(self.topoA.boundary), 9)
+        self.assertEqual(len(self.topoB.boundary), 6)
+        self.assertEqual(len(self.topoA.boundary['trimmed']), 3)
+        self.assertEqual(len(self.topoB.boundary['trimmed']), 3)
 
     def test_interfaces(self):
         self.assertEqual(len(self.topoA.interfaces), 4)
diff --git a/tests/test_points.py b/tests/test_points.py
index 80b53ee8b..30bec57e7 100644
--- a/tests/test_points.py
+++ b/tests/test_points.py
@@ -165,17 +165,16 @@ def test_weights(self):
             self.assertLess(abs(pnt.weights.sum()-exact), 1e-15)
 
     def test_points(self):
-        self.assertEqual(self.bezier.npoints, 27)
+        self.assertEqual(self.bezier.npoints, 26)
         for x in [0., .25, .5, .75, 1.]:
             for y in [0., .25, .5, .75, 1.]:
                 if x or y:
                     self.assertIn([x, y], self.bezier.coords.tolist())
         self.assertIn([0., .125], self.bezier.coords.tolist())
-        self.assertIn([.0625, .0625], self.bezier.coords.tolist())
         self.assertIn([.125, 0.], self.bezier.coords.tolist())
 
     def test_tri(self):
-        self.assertEqual(len(self.bezier.tri), 34)
+        self.assertEqual(len(self.bezier.tri), 33)
 
     def test_hull(self):
-        self.assertEqual(len(self.bezier.hull), 18)
+        self.assertEqual(len(self.bezier.hull), 17)

From bb992c335f7e4a8993da9c61ea32e7cc164fbd34 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Mon, 4 Jul 2022 16:41:27 +0200
Subject: [PATCH 18/19] remove cone, ribbons, edge_vertices base implem

This patch removes the Reference object's cone and ribbons attributes and
removes the naive base implementation of edge_vertices which is provided by all
the relevant remaining subclasses.
---
 nutils/element.py     | 148 ++----------------------------------------
 tests/test_element.py |   4 --
 tests/test_points.py  |  38 -----------
 3 files changed, 4 insertions(+), 186 deletions(-)

diff --git a/nutils/element.py b/nutils/element.py
index 61538167e..9d5655bdc 100644
--- a/nutils/element.py
+++ b/nutils/element.py
@@ -24,7 +24,7 @@ class Reference(types.Singleton):
     'reference element'
 
     __slots__ = 'ndims',
-    __cache__ = 'connectivity', 'edgechildren', 'ribbons', 'volume', 'centroid', '_linear_bernstein', 'getpoints'
+    __cache__ = 'connectivity', 'edgechildren', 'volume', 'centroid', '_linear_bernstein', 'getpoints'
 
     @types.apply_annotations
     def __init__(self, ndims: int):
@@ -83,14 +83,7 @@ def edge_vertices(self):
         `vol.vertices`).
         '''
 
-        # naive implementation; will be removed in a later commit
-        edge_vertices = []
-        for etrans, eref in self.edges:
-            dist = numpy.linalg.norm(self.vertices[:,_,:] - etrans.apply(eref.vertices), axis=2)
-            emap = numpy.argmin(dist, 0)
-            assert (dist[emap, numpy.arange(eref.nverts)] < 1e-15).all(), dist
-            edge_vertices.append(types.frozenarray(emap, copy=False))
-        return tuple(edge_vertices)
+        raise NotImplementedError
 
     __and__ = lambda self, other: self if self == other else NotImplemented
     __or__ = lambda self, other: self if self == other else NotImplemented
@@ -170,27 +163,6 @@ def edgechildren(self):
             edgechildren.append(types.frozenarray(children))
         return tuple(edgechildren)
 
-    @property
-    def ribbons(self):
-        # tuples of (iedge1,jedge1), (iedge2,jedge2) pairs
-        assert self.ndims >= 2
-        transforms = {}
-        ribbons = []
-        for iedge1, (etrans1, edge1) in enumerate(self.edges):
-            if edge1:
-                for iedge2, (etrans2, edge2) in enumerate(edge1.edges):
-                    if edge2:
-                        key = tuple(sorted(tuple(p) for p in (etrans1 * etrans2).apply(edge2.vertices)))
-                        try:
-                            jedge1, jedge2 = transforms.pop(key)
-                        except KeyError:
-                            transforms[key] = iedge1, iedge2
-                        else:
-                            assert self.edge_refs[jedge1].edge_refs[jedge2] == edge2
-                            ribbons.append(((iedge1, iedge2), (jedge1, jedge2)))
-        assert not transforms
-        return tuple(ribbons)
-
     def getischeme(self, ischeme):
         ischeme, degree = parse_legacy_ischeme(ischeme)
         points = self.getpoints(ischeme, degree)
@@ -290,9 +262,6 @@ def slice(self, levels, ndivisions):
 
         return MosaicReference(self, refs, midpoint)
 
-    def cone(self, trans, tip):
-        return Cone(self, trans, tip)
-
     def check_edges(self, tol=1e-15, print=print):
         volume = 0
         zero = 0
@@ -403,7 +372,7 @@ class SimplexReference(Reference):
     'simplex reference'
 
     __slots__ = ()
-    __cache__ = 'edge_refs', 'edge_transforms', 'edge_vertices', 'ribbons', '_get_poly_coeffs_bernstein', '_get_poly_coeffs_lagrange', '_integer_barycentric_coordinates'
+    __cache__ = 'edge_refs', 'edge_transforms', 'edge_vertices', '_get_poly_coeffs_bernstein', '_get_poly_coeffs_lagrange', '_integer_barycentric_coordinates'
 
     @property
     def vertices(self):
@@ -441,10 +410,6 @@ def child_refs(self):
     def child_transforms(self):
         return tuple(transform.SimplexChild(self.ndims, ichild) for ichild in range(2**self.ndims))
 
-    @property
-    def ribbons(self):
-        return tuple(((iedge1, iedge2), (iedge2+1, iedge1)) for iedge1 in range(self.ndims+1) for iedge2 in range(iedge1, self.ndims))
-
     def getpoints(self, ischeme, degree):
         if ischeme == 'gauss':
             return points.SimplexGaussPoints(self.ndims, degree if numeric.isint(degree) else sum(degree))
@@ -649,7 +614,7 @@ class TensorReference(Reference):
     'tensor reference'
 
     __slots__ = 'ref1', 'ref2'
-    __cache__ = 'vertices', 'edge_transforms', 'edge_vertices', 'ribbons', 'child_transforms', 'getpoints', 'get_poly_coeffs', 'centroid', 'simplices'
+    __cache__ = 'vertices', 'edge_transforms', 'edge_vertices', 'child_transforms', 'getpoints', 'get_poly_coeffs', 'centroid', 'simplices'
 
     def __init__(self, ref1, ref2):
         assert not isinstance(ref1, TensorReference)
@@ -756,29 +721,6 @@ def edge_refs(self):
             edge_refs.extend(self.ref1 * edge2 for edge2 in self.ref2.edge_refs)
         return tuple(edge_refs)
 
-    @property
-    def ribbons(self):
-        if self.ref1.ndims == 0:
-            return self.ref2.ribbons
-        if self.ref2.ndims == 0:
-            return self.ref1.ribbons
-        ribbons = []
-        for iedge1 in range(self.ref1.nedges):
-            #iedge = self.ref1.edge_refs[iedge] * self.ref2
-            for iedge2 in range(self.ref2.nedges):
-                #jedge = self.ref1 * self.ref2.edge_refs[jedge]
-                jedge1 = self.ref1.nedges + iedge2
-                jedge2 = iedge1
-                if self.ref1.ndims > 1:
-                    iedge2 += self.ref1.edge_refs[iedge1].nedges
-                ribbons.append(((iedge1, iedge2), (jedge1, jedge2)))
-        if self.ref1.ndims >= 2:
-            ribbons.extend(self.ref1.ribbons)
-        if self.ref2.ndims >= 2:
-            ribbons.extend(((iedge1+self.ref1.nedges, iedge2+self.ref1.nedges),
-                            (jedge1+self.ref1.nedges, jedge2+self.ref1.nedges)) for (iedge1, iedge2), (jedge1, jedge2) in self.ref2.ribbons)
-        return tuple(ribbons)
-
     @property
     def child_transforms(self):
         return tuple(transform.TensorChild(trans1, trans2) for trans1 in self.ref1.child_transforms for trans2 in self.ref2.child_transforms)
@@ -817,88 +759,6 @@ def _flat_refs(self):
                 yield ref
 
 
-class Cone(Reference):
-    'cone'
-
-    __slots__ = 'edgeref', 'etrans', 'tip', 'extnorm', 'height'
-    __cache__ = 'vertices', 'edge_transforms', 'edge_refs'
-
-    @types.apply_annotations
-    def __init__(self, edgeref, etrans, tip: types.arraydata):
-        assert etrans.fromdims == edgeref.ndims
-        assert etrans.todims == tip.shape[0]
-        super().__init__(etrans.todims)
-        self.edgeref = edgeref
-        self.etrans = etrans
-        self.tip = numpy.asarray(tip)
-        ext = etrans.ext
-        self.extnorm = numpy.linalg.norm(ext)
-        self.height = numpy.dot(etrans.offset - self.tip, ext) / self.extnorm
-        assert self.height >= 0, 'tip is positioned at the negative side of edge'
-
-    @property
-    def vertices(self):
-        return types.frozenarray(numpy.vstack([[self.tip], self.etrans.apply(self.edgeref.vertices)]), copy=False)
-
-    @property
-    def simplices(self):
-        return types.frozenarray([[0, *emap[simplex]] for emap, eref in zip(self.edge_vertices, self.edge_refs) for simplex in eref.simplices])
-
-    @property
-    def edge_transforms(self):
-        edge_transforms = [self.etrans]
-        if self.edgeref.ndims > 0:
-            for trans, edge in self.edgeref.edges:
-                if edge:
-                    b = self.etrans.apply(trans.offset)
-                    A = numpy.hstack([numpy.dot(self.etrans.linear, trans.linear), (self.tip-b)[:, _]])
-                    newtrans = transform.Updim(A, b, isflipped=self.etrans.isflipped ^ trans.isflipped ^ (self.ndims % 2 == 1))  # isflipped logic tested up to 3D
-                    edge_transforms.append(newtrans)
-        else:
-            edge_transforms.append(transform.Updim(numpy.zeros((1, 0)), self.tip, isflipped=not self.etrans.isflipped))
-        return tuple(edge_transforms)
-
-    @property
-    def edge_refs(self):
-        edge_refs = [self.edgeref]
-        if self.edgeref.ndims > 0:
-            extrudetrans = transform.Updim(numpy.eye(self.ndims-1)[:, :-1], numpy.zeros(self.ndims-1), isflipped=self.ndims % 2 == 0)
-            tip = numpy.array([0]*(self.ndims-2)+[1], dtype=float)
-            edge_refs.extend(edge.cone(extrudetrans, tip) for edge in self.edgeref.edge_refs if edge)
-        else:
-            edge_refs.append(getsimplex(0))
-        return tuple(edge_refs)
-
-    def getpoints(self, ischeme, degree):
-        if ischeme == 'gauss':
-            if self.nverts == self.ndims+1:  # use optimal gauss schemes for simplex-like cones
-                pnts = getsimplex(self.ndims).getpoints(ischeme, degree)
-                linearT = self.etrans.apply(self.edgeref.vertices) - self.tip
-                coords = pnts.coords @ linearT + self.tip
-                weights = pnts.weights * abs(numpy.linalg.det(linearT))
-                return points.CoordsWeightsPoints(coords, weights)
-            epoints = self.edgeref.getpoints('gauss', degree)
-            tx, tw = points.gauss((degree + self.ndims - 1)//2)
-            wx = tx**(self.ndims-1) * tw * self.extnorm * self.height
-            return points.CoordsWeightsPoints((tx[:, _, _] * (self.etrans.apply(epoints.coords)-self.tip)[_, :, :] + self.tip).reshape(-1, self.ndims), (epoints.weights[_, :] * wx[:, _]).ravel())
-        if ischeme in ('_centroid', 'uniform'):
-            layers = [(i+1, (i+.5)/degree) for i in range(degree)] if ischeme == 'uniform' \
-                else [(None, self.ndims/(self.ndims+1))]  # centroid
-            coords = numpy.concatenate([self.etrans.apply(self.edgeref.getpoints(ischeme, p).coords) * w + self.tip * (1-w) for p, w in layers])
-            volume = self.edgeref.volume * self.extnorm * self.height / self.ndims
-            return points.CoordsUniformPoints(coords, volume)
-        if ischeme == 'vtk' and self.nverts == 5 and self.ndims == 3:  # pyramid
-            return points.CoordsPoints(self.vertices[[1, 2, 4, 3, 0]])
-        return points.ConePoints(self.edgeref.getpoints(ischeme, degree), self.etrans, self.tip)
-
-    def inside(self, point, eps=0):
-        # point = etrans.apply(epoint) * xi + tip * (1-xi) => etrans.apply(epoint) = tip + (point-tip) / xi
-        xi = numpy.dot(self.etrans.ext, point-self.tip) / numpy.dot(self.etrans.ext, self.etrans.offset-self.tip)
-        return -eps <= xi <= 1+eps and self.edgeref.inside(numpy.linalg.solve(
-            numpy.dot(self.etrans.linear.T, self.etrans.linear),
-            numpy.dot(self.etrans.linear.T, self.tip + (point-self.tip)/xi - self.etrans.offset)), eps=eps)
-
-
 class OwnChildReference(Reference):
     'forward self as child'
 
diff --git a/tests/test_element.py b/tests/test_element.py
index f389dd250..5a599cfba 100644
--- a/tests/test_element.py
+++ b/tests/test_element.py
@@ -44,10 +44,6 @@ def test_swap(self):
                 swapped_down = etrans_.swapdown(ctrans_)
                 self.assertEqual(swapped_down, (etrans, ctrans))
 
-    @parametrize.enable_if(lambda ref, **kwargs: ref.ndims >= 2)
-    def test_ribbons(self):
-        self.ref.ribbons
-
     @parametrize.enable_if(lambda ref, **kwargs: not isinstance(ref, element.MosaicReference) and ref.ndims >= 1)
     def test_connectivity(self):
         for ichild, edges in enumerate(self.ref.connectivity):
diff --git a/tests/test_points.py b/tests/test_points.py
index 30bec57e7..3b215b07c 100644
--- a/tests/test_points.py
+++ b/tests/test_points.py
@@ -103,44 +103,6 @@ def test_pyramid(self):
                 self.assertIn(sorted(h), fullhull)
 
 
-@parametrize
-class cone(TestCase):
-
-    def setUp(self):
-        super().setUp()
-        if self.shape == 'square':
-            self.edgeref = element.getsimplex(1)**2
-        elif self.shape == 'triangle':
-            self.edgeref = element.getsimplex(2)
-        else:
-            raise Exception('invalid shape: {!r}'.format(self.shape))
-        self.etrans = transform.Updim(linear=[[-1., 0], [0, -3], [0, 0]], offset=[1., 3, 1], isflipped=False)
-        self.cone = element.Cone(edgeref=self.edgeref, etrans=self.etrans, tip=[1., 3, 0])
-
-    def test_volume(self):
-        numpy.testing.assert_almost_equal(actual=self.cone.volume, desired=self.edgeref.volume)
-
-    def _test_points(self, *args):
-        points = self.cone.getpoints(*args)
-        if hasattr(points, 'weights'):
-            numpy.testing.assert_almost_equal(actual=self.cone.volume, desired=points.weights.sum())
-        # check that all points lie within pyramid/prism
-        x, y, z = points.coords.T
-        self.assertTrue(numpy.all(numpy.greater_equal(x, 1-z) & numpy.less_equal(x, 1) & numpy.greater_equal(y, 1-z) & numpy.less_equal(y, 3)))
-        if self.shape == 'triangle':
-            self.assertTrue(numpy.less_equal(2-x-y/3, z).all())
-
-    def test_gauss(self):
-        self._test_points('gauss', 3)
-
-    def test_uniform(self):
-        self._test_points('uniform', 3)
-
-
-cone(shape='square')
-cone(shape='triangle')
-
-
 class trimmed(TestCase):
 
     def setUp(self):

From 232cad30d282afc67140a38b5919349fa3844433 Mon Sep 17 00:00:00 2001
From: Gertjan van Zwieten <gertjanvanzwieten@evalf.nl>
Date: Fri, 5 Aug 2022 16:31:06 +0200
Subject: [PATCH 19/19] add changelog entry

---
 CHANGELOG | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/CHANGELOG b/CHANGELOG
index 0ba43e533..436d5208d 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -5,6 +5,17 @@ The following overview lists user facing changes as well as newly added
 features in inverse chronological order.
 
 
+IMPROVED: more efficient trimming
+
+The trim routine (which is used for the Finite Cell Method) is rewritten for
+speed and to produce more efficient quadrature schemes. The changes relate to
+the subdivision at the deepest refinement level. While this step used to
+introduce auxiliary vertices at every dimension (lines, faces, volumes), the
+new implementation limits the introduction of vertices to the line segments
+only, resulting in a subdivision that consists of fewer simplices and
+consequently fewer quadrature points.
+
+
 REMOVED: Nutils configuration file
 
 Support for the Nutils configuration file (which used to be located in either