extend unit meshes (#798)

This PR extends `mesh.unitsquare` with a multipatch variant and adds
`mesh.unitcircle`.

examples/cahnhilliard.py

@@ -144,10 +144,11 @@ def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
 
     log.info('contact angle: {:.0f}°'.format(numpy.arccos((wtensn - wtensp) / stens) * 180 / numpy.pi))
 
-    domain, geom = mesh.unitsquare(nelems, etype)
     if circle:
-        angle = (geom-.5) * (numpy.pi/2)
-        geom = .5 + numpy.sin(angle) * numpy.cos(angle)[[1, 0]] / numpy.sqrt(2)
+        domain, geom = mesh.unitcircle(nelems, etype)
+        geom = (geom + 1) / 2
+    else:
+        domain, geom = mesh.unitsquare(nelems, etype)
 
     bezier = domain.sample('bezier', 5)  # sample for surface plots
     grid = domain.locate(geom, numeric.simplex_grid([1, 1], 1/40), maxdist=1/nelems, skip_missing=True, tol=1e-5)  # sample for quivers
@@ -242,16 +243,26 @@ def test_square(self):
                 eNoBYgCd/3TIkccBNkQ6IDqIN3/MF8cSx9Y02TmdOVHLMcecxxLIEjQUOAHOa8a1xWw3izb2M9UzPMc0
                 xmnGpzibODY34tETyJHHp8hbyWU2xzZTydfIOsrNyo3Gi8jCyyXIm8hkzD3K1IAxtQ==''')
 
-    def test_mixedcircle(self):
+    def test_multipatchcircle(self):
         args = main(size=parse('10cm'), epsilon=parse('5cm'), mobility=parse('1μL*s/kg'),
                      stens=parse('50mN/m'), wtensn=parse('30mN/m'), wtensp=parse('20mN/m'), nelems=3,
-                     etype='mixed', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
+                     etype='multipatch', degree=2, timestep=parse('1h'), tol=parse('1nJ/m'),
                      endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True)
         with self.subTest('concentration'):
             self.assertAlmostEqual64(args['φ'], '''
-                eNoBYgCd/w01AjX+NAw1IjXTNMw0iTRPNDI1vDQcNTk0uzJ9NFM0HS4P0SbMcssOy0wzZjNw0b0zljHK
-                z6U0ps8zM/LPjspVypDKUsuLzk3MgM3OzYnN7s/61KfP2zH4MADNhst3z7DMoBcvyQ==''')
+                eNoNz01IVGEUxnFBSgxbVIqCCEJQi4E7933PcVGL/FrISNDGFsUo5ORCBovQIctxQBQRZFpFIvZJ4aYQ
+                YQQlkgTL1XPOe0dkKEOhhSJkOz8KRL38F8/2+Tl2/JW/hL3icj7PSjOc4CdczMO0ZEv5j22nKZMy6VVD
+                Z+mh3TBp0+hnzQOvJZJwA65E/0mbVK89zovbdftu2G3qjlx1Ma2XX3gjZ2RWcvJMshhFH1rkoqwjKVek
+                E234gAy2sIIY2hFHB15qnb6XEfPUDPqeOPheyg2hC5BDfKQc9dFB2DpdpoJZpCr6++On/DcNwRzP8SQf
+                0R6tUZprOEdN9I7q7LTt9itoO9xFc9sk/Q47bV544/mkXwgS7oYe5xuCGhfX73JXItLqTjSj5R4XeoIK
+                HZFKV6QXZB4xUfSaInM9eo6O7JgtW13WXdMcfS47iLt+OQw/faZH1GlfR2/p22hz8Fu/SZVGpEw+IRu6
+                amUBC6H/Hm7iDlK4hgn04H5onsElOQVf7chD''')
         with self.subTest('chemical-potential'):
             self.assertAlmostEqual64(args['η'], '''
-                eNoBYgCd/+s1Bcp4ztI3gjYFyZk4YzVjyfA2AzdAMj032zfLNTE4fMm7yLnGisbqxZPJ2MsfyD81csiv
-                x+E5xDhjOJA3msZ1xZTFa8ddx/fG88eCx73H1MieM/c0WDihMUrLvMYZNpvIrWQ0sw==''')
+                eNoN0MtLlGEUB2AjoiIIKzChXS5cVOR837znHMNAGaaLQWAxi7YudOFl4WiWRchsXIhQXsAYKlIIjGiS
+                TDeRs6nfOe/7jUUWuWnjostGwwoXbZLnP3iykhUnW7xPqmSOlilLdbJLfvMkz7g49ZfreS/fTT/3jfbK
+                zbmN1Kr74m5Gv/wP/ehf6k+UsIg8Lsfr0eMEvuIT3+4f6qofthIuYVyXcdjX+gP+qnZqi6qVbE23bI8d
+                1Hk06GtktFp7cQ+9KOBIshEWrY+e0gma0X/od0vJGIYwYCm9Ly+kKCs7ylKQJV6QjFTRhD3h26ebuYcf
+                cDG5UplNf+b90kUnuUnaeZ7H6aI/7t+GGtfmvsb5KOcuvF8P3WTpY+GOz9Ah6khfS27ZiE5qMYyEWT8d
+                jUb5yrZ907PhupE+QkHfoSZei3hlgT5Rgd4Eb6eoHJf1O46GCW2VTmkTkn7+4z5Yl0tiZ9uaS6ptWnfr
+                JqYwiBs4gyw6kMM59KF1Z60WDWjGM5zX/wCZyEY=''')

nutils/mesh.py

@@ -629,6 +629,19 @@ def unitsquare(nelems, etype):
 
         * ``"mixed"``: unstructured mesh of triangles and squares.
 
+        * ``"multipatch"``: multipatch mesh consisting of five patches with the
+            following structure:
+
+            .. code-block:: none
+
+               *─────*
+               │╲   ╱│
+               │ *─* │
+               │ │ │ │
+               │ *─* │
+               │╱   ╲│
+               *─────*
+
     Returns
     -------
     :class:`nutils.topology.TransformChainsTopology`:
@@ -641,6 +654,7 @@ def unitsquare(nelems, etype):
 
     if etype == 'square':
         topo, geom = rectilinear([nelems, nelems], space=space)
+        return topo, geom / nelems
 
     elif etype in ('triangle', 'mixed'):
         simplices = numpy.concatenate([
@@ -672,11 +686,72 @@ def unitsquare(nelems, etype):
         x, y = topo.boundary.sample('_centroid', None).eval(geom).T
         btopo = topo.boundary
         topo = topo.withboundary(left=btopo[x<.1], right=btopo[x>nelems-.1], bottom=btopo[y<.1], top=btopo[y>nelems-.1])
+        return topo, geom / nelems
+
+    elif etype == 'multipatch':
+        # 2─────3
+        # │╲   ╱│
+        # │ 6─7 │
+        # │ │ │ │
+        # │ 4─5 │
+        # │╱   ╲│
+        # 0─────1
+        topo, geom = multipatch(
+            patches=[[0, 4, 1, 5], [2, 6, 3, 7], [0, 4, 2, 6], [1, 5, 3, 7], [4, 6, 5, 7]],
+            patchverts=[[0, 0], [3, 0], [0, 3], [3, 3], [1, 1], [2, 1], [1, 2], [2, 2]],
+            nelems=nelems)
+        topo = topo.withboundary(
+            bottom=topo['patch0'].boundary['bottom'],
+            top=topo['patch1'].boundary['bottom'],
+            left=topo['patch2'].boundary['bottom'],
+            right=topo['patch3'].boundary['bottom'])
+        return topo, geom / 3
 
     else:
         raise Exception('invalid element type {!r}'.format(etype))
 
-    return topo, geom / nelems
+
+def _square_to_circle(geom):
+    angle = (geom - 0.5) * (numpy.pi / 2)
+    return numpy.sqrt(2) * numpy.sin(angle) * numpy.cos(angle)[[1, 0]]
+
+
+def unitcircle(nelems: int, variant: str) -> Tuple[Topology, function.Array]:
+    '''Unit circle mesh.
+
+    The circle is centered at coordinate ``(0, 0)`` and has unit radius.
+
+    Args
+    ----
+    nelems : :class:`int`
+        Number of elements along boundary
+    variant : :class:`str`
+        Mesh variant. Supported are:
+
+        *   ``"rectilinear"``: structured mesh of squares, blown-up to a circle.
+
+        *   ``"multipatch"``: multipatch mesh consisting of five patches with
+            the same topological structure as :func:`unitsquare` with
+            ``etype='multipatch'``.
+
+    Returns
+    -------
+    :class:`nutils.topology.TransformChainsTopology`:
+        The structured/unstructured topology.
+    :class:`nutils.function.Array`:
+        The geometry function.
+    '''
+
+    if variant == 'rectilinear':
+        topo, geom = unitsquare(nelems, 'square')
+        return topo, _square_to_circle(geom)
+
+    elif variant == 'multipatch':
+        topo, geom = unitsquare(nelems, 'multipatch')
+        return topo, _square_to_circle(geom)
+
+    else:
+        raise Exception('invalid variant {!r}'.format(variant))
 
 
 # vim:sw=4:sts=4:et

tests/test_mesh.py

@@ -144,3 +144,59 @@ def test_pum(self):
                     basis = self.domain.basis(basistype, degree=degree)
                     values = self.domain.interfaces.sample('uniform', 2).eval(basis*function.J(self.geom))
                     numpy.testing.assert_almost_equal(values.sum(1), 1)
+
+
+@parametrize
+class unitsquare(TestCase):
+
+    def setUp(self):
+        super().setUp()
+        self.domain, self.geom = mesh.unitsquare(nelems=4, etype=self.etype)
+
+    def test_volume(self):
+        self.assertAllAlmostEqual(self.domain.volume(self.geom), 1)
+
+    def test_boundaries(self):
+        self.assertAllAlmostEqual(self.domain.boundary.volume(self.geom), 4)
+        self.domain.check_boundary(geometry=self.geom)
+
+    def test_boundary_groups(self):
+        numpy.testing.assert_almost_equal(self.domain.boundary['left'].sample('gauss', 0).eval(self.geom[0]), 0)
+        numpy.testing.assert_almost_equal(self.domain.boundary['bottom'].sample('gauss', 0).eval(self.geom[1]), 0)
+        numpy.testing.assert_almost_equal(self.domain.boundary['right'].sample('gauss', 0).eval(self.geom[0]), 1)
+        numpy.testing.assert_almost_equal(self.domain.boundary['top'].sample('gauss', 0).eval(self.geom[1]), 1)
+
+    def test_interface(self):
+        geomerr = self.domain.interfaces.sample('uniform', 2).eval(self.geom - function.opposite(self.geom))
+        numpy.testing.assert_almost_equal(geomerr, 0, decimal=15)
+        normalerr = self.domain.interfaces.sample('uniform', 2).eval(self.geom.normal() + function.opposite(self.geom.normal()))
+        numpy.testing.assert_almost_equal(normalerr, 0, decimal=15)
+
+unitsquare(etype='square')
+unitsquare(etype='triangle')
+unitsquare(etype='mixed')
+unitsquare(etype='multipatch')
+
+
+@parametrize
+class unitcircle(TestCase):
+
+    def setUp(self):
+        super().setUp()
+        self.domain, self.geom = mesh.unitcircle(nelems=8, variant=self.variant)
+
+    def test_volume(self):
+        self.assertAllAlmostEqual(self.domain.volume(self.geom, degree=6), numpy.pi)
+
+    def test_boundaries(self):
+        self.assertAllAlmostEqual(self.domain.boundary.volume(self.geom, degree=6), 2 * numpy.pi)
+        self.domain.check_boundary(geometry=self.geom, degree=8)
+
+    def test_interface(self):
+        geomerr = self.domain.interfaces.sample('uniform', 2).eval(self.geom - function.opposite(self.geom))
+        numpy.testing.assert_almost_equal(geomerr, 0, decimal=15)
+        normalerr = self.domain.interfaces.sample('uniform', 2).eval(self.geom.normal() + function.opposite(self.geom.normal()))
+        numpy.testing.assert_almost_equal(normalerr, 0, decimal=15)
+
+unitcircle(variant='rectilinear')
+unitcircle(variant='multipatch')