From a041fc662ad18f3b8c7ac4aa5f80ee0c45071d85 Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Wed, 20 Dec 2023 10:43:23 +0100 Subject: [PATCH 1/3] rewrite Topology.trim without evaluable.Points This patch rewrites `Topology.trim` with a `leveltopo` without `evaluable.Points`. Instead, the `levelset` is lowered with points sequence vertex for all relevant degrees. --- nutils/element.py | 6 +++--- nutils/topology.py | 24 ++++++++++++++++++------ 2 files changed, 21 insertions(+), 9 deletions(-) diff --git a/nutils/element.py b/nutils/element.py index 3f8787f9b..34ac3621e 100644 --- a/nutils/element.py +++ b/nutils/element.py @@ -282,7 +282,7 @@ def _linear_cover(self, ctransforms, maxrefine): if len(ctransforms) == 1: ctrans, = ctransforms assert not ctrans - return ((), self.getpoints('vertex', maxrefine), allindices), + return ((), maxrefine, allindices), if maxrefine == 0: raise Exception('maxrefine is too low') cbins = [set() for ichild in range(self.nchildren)] @@ -292,9 +292,9 @@ def _linear_cover(self, ctransforms, maxrefine): if not all(cbins): raise Exception('transformations to not form an element cover') fcache = cache.WrapperCache() - return tuple(((ctrans,) + trans, points, cindices[indices]) + return tuple(((ctrans,) + trans, degree, 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._linear_cover](frozenset(cbin), maxrefine-1)) + for trans, degree, 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 a28dff5df..0f2a7ad90 100644 --- a/nutils/topology.py +++ b/nutils/topology.py @@ -1471,10 +1471,18 @@ def trim(self, levelset, maxrefine, ndivisions=8, name='trimmed', leveltopo=None levels = levelset.eval(_trim_index=ielem, **arguments) refs.append(ref.trim(levels, maxrefine=maxrefine, ndivisions=ndivisions)) else: + # `levelset` is evaluable on `leveltopo`, which must be a uniform + # or hierarchical refinement of `self`. Each element of `self` is + # trimmed given a sample of `levelset` at the pointset `vertex` + # with degree `maxrefine`. This sample is obtained by evaluating + # `levelset` for each child element in `leveltopo` at the relevant + # subset of the parent sample. Per child element this subset is the + # pointset `vertex` with degree `maxrefine` minus the ancestral + # level. The mapping from child points to parent points and the + # degree of the pointset is provided by `Element._linear_cover`. log.info('collecting leveltopo elements') - coordinates = evaluable.Points(evaluable.NPoints(), evaluable.constant(self.ndims)) - ielem = evaluable.Argument('_leveltopo_ielem', (), int) - levelset = levelset.lower(function.LowerArgs.for_space(self.space, (leveltopo.transforms, leveltopo.opposites), ielem, coordinates)).optimized_for_numpy + lowered_levelset = {} + ielem_arg = evaluable.InRange(evaluable.Argument('_ielem', (), int), evaluable.constant(len(leveltopo))) bins = [set() for ielem in range(len(self))] for trans in leveltopo.transforms: ielem, tail = self.transforms.index_with_tail(trans) @@ -1487,10 +1495,14 @@ def trim(self, levelset, maxrefine, ndivisions=8, name='trimmed', leveltopo=None # confirm cover and greedily optimize order mask = numpy.ones(len(levels), dtype=bool) while mask.any(): - imax = numpy.argmax([mask[indices].sum() for tail, points, indices in cover]) - tail, points, indices = cover.pop(imax) + imax = numpy.argmax([mask[indices].sum() for tail, degree, indices in cover]) + tail, degree, indices = cover.pop(imax) ielem = leveltopo.transforms.index(trans + tail) - levels[indices] = levelset.eval(_leveltopo_ielem=ielem, _points=points, **arguments) + if degree not in lowered_levelset: + coordinates = leveltopo.references.getpoints('vertex', degree).get_evaluable_coords(ielem_arg) + lower_args = function.LowerArgs.for_space(self.space, (leveltopo.transforms, leveltopo.opposites), ielem_arg, coordinates) + lowered_levelset[degree] = levelset.lower(lower_args).optimized_for_numpy + levels[indices] = lowered_levelset[degree].eval(_ielem=ielem, **arguments) mask[indices] = False refs.append(ref.trim(levels, maxrefine=maxrefine, ndivisions=ndivisions)) log.debug('cache', fcache.stats) From 477d006a9494775585b21ebc08b589620b4683ef Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Wed, 20 Dec 2023 11:00:48 +0100 Subject: [PATCH 2/3] rewrite unittest without evaluable.Points --- tests/test_topology.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/tests/test_topology.py b/tests/test_topology.py index ab4e0920a..2fe6aaebb 100644 --- a/tests/test_topology.py +++ b/tests/test_topology.py @@ -649,9 +649,9 @@ def assertConnectivity(self, domain, geom): interfaces = domain.interfaces bmask = numpy.zeros(len(boundary), dtype=int) imask = numpy.zeros(len(interfaces), dtype=int) - coordinates = evaluable.Points(evaluable.NPoints(), evaluable.constant(boundary.ndims)) edges = domain.transforms.edges(domain.references) iedge = evaluable.Argument('_iedge', (), int) + coordinates = domain.references.edges.getpoints('gauss', 2).get_evaluable_coords(iedge) lowered_geom = geom.lower(function.LowerArgs.for_space(domain.space, (edges,), iedge, coordinates)).simplified for ielem, ioppelems in enumerate(domain.connectivity): for iedge, ioppelem in enumerate(ioppelems): @@ -673,9 +673,8 @@ def assertConnectivity(self, domain, geom): self.assertEqual(interfaces.opposites[index], opptrans) imask[index] += 1 self.assertEqual(eref, opperef) - points = eref.getpoints('gauss', 2) - a0 = lowered_geom.eval(_iedge=edges.index(trans), _points=points) - a1 = lowered_geom.eval(_iedge=edges.index(opptrans), _points=points) + a0 = lowered_geom.eval(_iedge=edges.index(trans)) + a1 = lowered_geom.eval(_iedge=edges.index(opptrans)) numpy.testing.assert_array_almost_equal(a0, a1) self.assertTrue(numpy.equal(bmask, 1).all()) self.assertTrue(numpy.equal(imask, 2).all()) From 3aa7a6fa4ec8be099af237ad4fe5a9da9dff8c01 Mon Sep 17 00:00:00 2001 From: Joost van Zwieten Date: Wed, 20 Dec 2023 11:06:47 +0100 Subject: [PATCH 3/3] remove unused evaluable.Points,NPoints,Weights --- nutils/evaluable.py | 37 ------------------------------------- tests/test_evaluable.py | 3 --- 2 files changed, 40 deletions(-) diff --git a/nutils/evaluable.py b/nutils/evaluable.py index 1b91ee7f6..a138032d0 100644 --- a/nutils/evaluable.py +++ b/nutils/evaluable.py @@ -1024,43 +1024,6 @@ def _const_uniform(self): return lower if lower == upper else None -class NPoints(Array): - 'The length of the points axis.' - - def __init__(self): - super().__init__(args=(EVALARGS,), shape=(), dtype=int) - - @staticmethod - def evalf(evalargs): - points = evalargs['_points'].coords - return types.frozenarray(points.shape[0]) - - def _intbounds_impl(self): - return 0, float('inf') - - -class Points(Array): - - def __init__(self, npoints, ndim): - super().__init__(args=(EVALARGS,), shape=(npoints, ndim), dtype=float) - - @staticmethod - def evalf(evalargs): - return evalargs['_points'].coords - - -class Weights(Array): - - def __init__(self, npoints): - super().__init__(args=(EVALARGS,), shape=(npoints,), dtype=float) - - @staticmethod - def evalf(evalargs): - weights = evalargs['_points'].weights - assert numeric.isarray(weights) and weights.ndim == 1 - return weights - - class Orthonormal(Array): 'make a vector orthonormal to a subspace' diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py index cadab5b6f..1da113e27 100644 --- a/tests/test_evaluable.py +++ b/tests/test_evaluable.py @@ -755,9 +755,6 @@ def test_inrange_strict(self): def test_inrange_empty(self): self.assertEqual(evaluable.InRange(self.S('n', float('-inf'), float('inf')), evaluable.constant(0))._intbounds, (0, 0)) - def test_npoints(self): - self.assertEqual(evaluable.NPoints()._intbounds, (0, float('inf'))) - def test_bool_to_int(self): self.assertEqual(evaluable.BoolToInt(evaluable.constant(numpy.array([False, True], dtype=bool)))._intbounds, (0, 1))