Fix "not watertight" error in 3D trimming operation

nutils/element.py

@@ -185,8 +185,7 @@
 
     @cached_property
     def volume(self):
-        volume, = self.getpoints('_centroid', None).weights
-        return volume
+        return self.getpoints('gauss', 1).weights.sum()
 
     @cached_property
     def centroid(self):
@@ -211,17 +210,14 @@
         # slice along levelset by recursing over dimensions
 
         assert len(levels) == len(self.vertices)
-        if numpy.greater_equal(levels, 0).all():
+        L = levels[self.simplices]
+        if numpy.greater_equal(L, 0).all():
             return self
-        if numpy.less_equal(levels, 0).all():
+        if numpy.less_equal(L, 0).all():
             return self.empty
         assert self.ndims >= 1
 
         refs = tuple(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
 
         if self.ndims == 1:
 
@@ -232,6 +228,7 @@
 
             iedge = [i for (i,), edge in zip(self.edge_vertices, self.edge_refs) if edge]
             l0, l1 = levels[iedge]
+            assert l0 * l1 < 0
             nbins = 2**ndivisions
             xi = numpy.round(l0/(l0-l1) * nbins)
             if xi in (0, nbins):
@@ -244,20 +241,33 @@
             # 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])
+            # fully retained or fully removed) then we first consider if the
+            # hull amounts to the full or empty element; if not 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, newedge in zip(self.edge_transforms, self.edge_refs, refs):
+                if newedge != edge and isinstance(newedge, MosaicReference):
+                    midpoint = trans.apply(numpy.asarray(newedge._midpoint))
                     break
-            else:
+            else: # no new edge intersections
+                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
                 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]
 
+            # The reason that the two returns above are restricted to the
+            # situation that the levelset values did not cause any slicing
+            # between vertices is that these intersections might continue on an
+            # adjoining edge, so that local interventions lead to _ribs
+            # conflicts. The consequence of this is that the resulting mosaic
+            # may have a volume that is equal to zero or self.
+
         return MosaicReference(self, refs, types.arraydata(midpoint))
 
     def check_edges(self, tol=1e-15, print=print):
@@ -310,6 +320,35 @@
     def get_edge_dofs(self, degree, iedge):
         raise NotImplementedError
 
+    @cached_property
+    def _ribs(self):
+        # Index pairs of the form (i1,j1), (i2,j2), such that edge j1 of edge
+        # i1 coincides for edge j2 of edge i1. To identify matching edge-edges
+        # we map their vertices to the numbering of the baseref for comparison.
+        # Since by construction all coinciding edge-edges have equal references
+        # and orientation, we can do the identification on the basis of
+        # edge_vertices alone. We include the reference objects for good
+        # measure, to ensure that the conventions are upheld.
+        #
+        # NOTE: This attribute is used only by MosaicReference.__init__, but it
+        # is defined here so that it can be cached with the element.
+
+        pairs = []
+        seen = {}
+        for i, (erefi, emapi) in enumerate(zip(self.edge_refs, self.edge_vertices)):
+            if not isinstance(erefi, EmptyLike):
+                for j, (erefj, emapj) in enumerate(zip(erefi.edge_refs, erefi.edge_vertices)):
+                    if not isinstance(erefj, EmptyLike):
+                        key = erefj, tuple(emapi[emapj])
+                        try:
+                            i2, j2 = seen.pop(key)
+                        except KeyError:
+                            seen[key] = i, j
+                        else: # a counterpart is found, placing newedge2 against newedge2_
+                            pairs.append(((i, j), (i2, j2)))
+        assert not seen, f'leftover unmatched edges ({seen}) indicate the edges of {self} are not watertight!'
+        return tuple(pairs)
+
 
 class EmptyLike(Reference):
     'inverse reference element'
@@ -893,16 +932,20 @@
         # edges of the mosaic, and setting up the corresponding edge-vertex
         # relations.
 
+        edge_vertices = []
+        orientation = []
+
         if baseref.ndims == 1:
 
             # 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.
 
             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'
+            edge_vertices.extend(baseref.edge_vertices)
+            edge_vertices.append(types.frozenarray([imidpoint]))
+            isflipped, = [not trans.isflipped for trans, edge in zip(baseref.edge_transforms, edge_refs) if edge]
+            orientation.append(isflipped)
 
         else:
 
@@ -913,7 +956,6 @@
             # 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:]):
                     for i, v_ in enumerate(vertices[baseref.nverts:], start=baseref.nverts):
@@ -931,8 +973,11 @@
             # 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))
+            trimmed = []
+            for i, (edge, newedge) in enumerate(zip(baseref.edge_refs, edge_refs)):
+                n = edge.nedges
+                assert edge.edge_transforms == newedge.edge_transforms[:n]
+                trimmed.extend((i, j) for j in range(n, newedge.nedges))
 
             # The latter, however, is more tricky. This is the situation that
             # occurs, for instance, when two out of four edges of a square are
@@ -941,66 +986,32 @@
             # 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
+            # counterparts (cached in the _ribs attribute of the base
+            # reference) 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]]:
-                        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!'
+            for (i1, j1), (i2, j2) in baseref._ribs:
+                e1 = edge_refs[i1].edge_refs[j1]
+                e2 = edge_refs[i2].edge_refs[j2]
+                if e1 and not e2:
+                    trimmed.append((i1, j1))
+                elif e2 and not e1:
+                    trimmed.append((i2, j2))
+                elif e1 != e2:
+                    raise NotImplementedError
+
+            # What remains is only to extend the edge-vertex relations and
+            # to track if the new edges are left- or right-handed.
+
+            for i, j in trimmed:
+                newedge = edge_refs[i]
+                if not isinstance(newedge.edge_refs[j], SimplexReference):
+                    raise NotImplementedError
+                emap = edge_vertices[i][newedge.edge_vertices[j]]
+                if imidpoint not in emap:
+                    edge_vertices.append(types.frozenarray([imidpoint, *emap]))
+                    orientation.append(not baseref.edge_transforms[i].isflipped ^ newedge.edge_transforms[j].isflipped)
 
         self.baseref = baseref
         self._edge_refs = edge_refs
@@ -1063,6 +1074,8 @@
     def getpoints(self, ischeme, degree):
         if ischeme == 'vertex':
             return self.baseref.getpoints(ischeme, degree)
+        elif len(self.simplices) == 0:
+            return points.CoordsWeightsPoints(types.arraydata(numpy.empty((0,self.ndims))), types.arraydata(numpy.empty((0,))))
         elif ischeme in ('gauss', 'uniform', 'bezier'):
             simplexpoints = getsimplex(self.ndims).getpoints(ischeme, degree)
             subpoints = tuple(points.TransformPoints(simplexpoints, strans) for strans in self.simplex_transforms)

tests/test_finitecell.py

@@ -1,4 +1,4 @@
-from nutils import topology, mesh, function, evaluable
+from nutils import topology, mesh, function, evaluable, element
 from nutils.testing import TestCase, parametrize
 import treelog as log
 import numpy
@@ -159,6 +159,18 @@ def test_inter_elem(self):
                 with self.subTest(how=how, maxrefine=maxrefine):
                     self.domain.trim(z-.75+self.eps*perturb, maxrefine=maxrefine).check_boundary(self.geom, elemwise=True, print=self.fail)
 
+    def test_rib_mosaic(self):
+        baseref = element.getsimplex(3)
+        pos = baseref.trim([-7.68105356e-02, -4.28304119e-02, 9.65282734e-05, 1.57188628e-02], maxrefine=0, ndivisions=8)
+        # These weights cause one of the edges to vanish entirely, while a
+        # neighbouring edge is trimmed at the joint rib. The incompatibility is
+        # resolved by retaining the mosaic as a degenerate element.
+        assert any(isinstance(edge, element.MosaicReference) and edge.volume == 0 for edge in pos.edge_refs)
+        neg = baseref - pos
+        for ref in pos, neg:
+            ref._ribs # check consistency of the new mosaic
+            ref.check_edges(print=self.fail)
+
 
 class setoperations(TestCase):