From adcd5d7963ecea6b57bc1308f048e89e6e37163c Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Tue, 22 Aug 2023 11:19:23 +0200 Subject: [PATCH 1/7] fix faulty assertion: nverts should be nedges This patch fixes an assertion in the MosaicReference constructor that checks that the edge transforms of the new edges coinside with the edge transforms of the edges of the base reference. Here `nedges` was mistaken for `nverts`, which worked because in dimensions 1 and 2 the two are often equal, but which was technically incorrect. --- nutils/element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nutils/element.py b/nutils/element.py index 34ac3621e..4f588c714 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -931,7 +931,7 @@ def __init__(self, baseref: Reference, edge_refs: Tuple[Reference,...], midpoint # 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] + assert all(edge.edge_transforms == newedge.edge_transforms[:edge.nedges] for edge, newedge in zip(baseref.edge_refs, edge_refs)) # The latter, however, is more tricky. This is the situation that From ebd6c213304fa858ac0a176604882ad13d122f6c Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 23 Aug 2023 16:13:37 +0200 Subject: [PATCH 2/7] compute element volume based on gaus1 i/o centroid This patch changes the implementation of Element.volume to use degree=1 gauss quadrature points, rather than _centroid. Though the two are often equal, the new implementation is robust against the situation that an element has zero quadrature points. It also avoids the unnecessary computation of the centroid coordinate which is not used for the volume. --- nutils/element.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/nutils/element.py b/nutils/element.py index 4f588c714..aab978e3b 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -185,8 +185,7 @@ def with_children(self, child_refs): @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): From 26ac3b747959434c8ce7b08cd37b5fd4897e59a2 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 23 Aug 2023 16:16:20 +0200 Subject: [PATCH 3/7] fix getpoints for empty mosaic elements This patch adds a special case to MosaicReference.getpoints for the situation that its simplices are empty, in which case ConcatPoints crashed as it tried to establish the dimension. The new implementation creates a CoordsWeightsPoints instance of zero points. --- nutils/element.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/nutils/element.py b/nutils/element.py index aab978e3b..7b53c2703 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -1062,6 +1062,8 @@ def simplices(self): 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) From 4f1a4334f8798a4e5fd65b099d4debf0e4339136 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 23 Aug 2023 17:01:45 +0200 Subject: [PATCH 4/7] improve early return in Reference.slice This patch changes the sign test in Reference.slice that results in an early return if the levelset values are all of the same sign. By focusing this tests on the vertices used by simplices, rather than all vertices (some of which may not play a role due to earlier slicing) more cases might be caught here rather than in the later edge test, which is now restricted to ndims >= 2. --- nutils/element.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/nutils/element.py b/nutils/element.py index 7b53c2703..2a53da48f 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -210,17 +210,14 @@ def slice(self, levels, ndivisions): # 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: @@ -231,6 +228,7 @@ def slice(self, levels, ndivisions): 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): @@ -240,6 +238,11 @@ def slice(self, levels, ndivisions): else: + 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 + # 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 From c96298a9c0e51862d9db55f3179a070e4fb1ce31 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 23 Aug 2023 17:08:41 +0200 Subject: [PATCH 5/7] fix 'not watertight' trimming issues This patch fixes a ndivisions-related trimming issue where, due to rouding, adjoining edges of a newly formed mosaic are incompatible. The patch fixes this by retaining degenerate mosaics if they serve the purpose of completing the rib skeleton. --- nutils/element.py | 34 +++++++++++++++++++++------------- tests/test_finitecell.py | 13 ++++++++++++- 2 files changed, 33 insertions(+), 14 deletions(-) diff --git a/nutils/element.py b/nutils/element.py index 2a53da48f..01ebc3520 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -238,28 +238,36 @@ def slice(self, levels, ndivisions): else: - 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 - # 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): diff --git a/tests/test_finitecell.py b/tests/test_finitecell.py index 4bf90f904..fd76d4b8d 100644 --- a/tests/test_finitecell.py +++ b/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,17 @@ 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.check_edges(print=self.fail) + class setoperations(TestCase): From 5ee6afcd55d896adb8680a58cdf2c7a3e7a485cd Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 21 Aug 2023 21:11:19 +0200 Subject: [PATCH 6/7] introduce Reference._ribs This patch introduces the cached Reference._ribs property that allows an element's edge-edge relations to be reused in the formation of new mosaic elements. --- nutils/element.py | 113 +++++++++++++++++++-------------------- tests/test_finitecell.py | 1 + 2 files changed, 57 insertions(+), 57 deletions(-) diff --git a/nutils/element.py b/nutils/element.py index 01ebc3520..b4b228579 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -320,6 +320,35 @@ def get_poly_coeffs(self, basis, **kwargs): 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' @@ -941,8 +970,11 @@ def __init__(self, baseref: Reference, edge_refs: Tuple[Reference,...], midpoint # trivial, following the convention that existing edge transforms # are copied over in the modified edge. - assert all(edge.edge_transforms == newedge.edge_transforms[:edge.nedges] - 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 @@ -951,66 +983,33 @@ def __init__(self, baseref: Reference, edge_refs: Tuple[Reference,...], midpoint # 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. + 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 - 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. + # What remains is only to extend the edge-vertex relations and + # to track if the new edges are left- or right-handed. - 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!' + orientation = [] + 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 diff --git a/tests/test_finitecell.py b/tests/test_finitecell.py index fd76d4b8d..aca292680 100644 --- a/tests/test_finitecell.py +++ b/tests/test_finitecell.py @@ -168,6 +168,7 @@ def test_rib_mosaic(self): 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) From d7e8ff8806ffcbfbd31ad4c040b08baa3685302c Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Wed, 23 Aug 2023 16:21:32 +0200 Subject: [PATCH 7/7] improve organisation of the mosaic constructor This patch declares the edge_vertices and orientation lists to outside the if-block in which they are populated, for reasons of readibility. --- nutils/element.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/nutils/element.py b/nutils/element.py index b4b228579..fcd2f6560 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -932,6 +932,9 @@ def __init__(self, baseref: Reference, edge_refs: Tuple[Reference,...], midpoint # 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 @@ -939,9 +942,10 @@ def __init__(self, baseref: Reference, edge_refs: Tuple[Reference,...], midpoint # 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: @@ -952,7 +956,6 @@ def __init__(self, baseref: Reference, edge_refs: Tuple[Reference,...], midpoint # 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): @@ -1001,7 +1004,6 @@ def __init__(self, baseref: Reference, edge_refs: Tuple[Reference,...], midpoint # What remains is only to extend the edge-vertex relations and # to track if the new edges are left- or right-handed. - orientation = [] for i, j in trimmed: newedge = edge_refs[i] if not isinstance(newedge.edge_refs[j], SimplexReference):