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)