Unitcircle multipatch (#802)

This PR replaces the sin/cos based square to circle mapping for the
multipatch variant of the unit circle by a NURBS based approach. The new
map has the benefit of having a strictly positive Jacobian throughout
the domain, including everywhere on the boundary.

examples/cahnhilliard.py

@@ -115,7 +115,7 @@ def main(size: Length, epsilon: Length, mobility: Time/Density, stens: Tension,
        nelems [0]
          Number of elements along domain edge. When set to zero a value is set
          automatically based on the configured domain size and epsilon.
-       etype [square]
+       etype [rectilinear]
          Type of elements (square/triangle/mixed).
        degree [2]
          Polynomial degree.
@@ -250,19 +250,19 @@ def test_multipatchcircle(self):
                      endtime=parse('2h'), seed=0, circle=True, stab=stab.linear, showflux=True)
         with self.subTest('concentration'):
             self.assertAlmostEqual64(args['φ'], '''
-                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''')
+                eNoNyD9IlHEYB3BIUyI0HaSkyEkl4nzf3/M8EDk4Zw5x4BK5RAbRICWIUMshoaKBEaZ3mprhUi29w3ku
+                ZUtH0PN93p9/eBvuriWEMghsuCk0x8+nIolsyD85kIKMyl05LXNyVR5KhXO8SZ+4SHU87brdyOWvvMq1
+                rDROOddC1eBNUPR9Pm8TtoeRcDW4mfy41LN7favZT1mbL1sBt9BgGZyId23Y+hGiBmwJbuOcvUKikR5p
+                UbtRjxea1Xmd1SE/EHfYGdpxmXAQVX2Z7Nva8S+gFRG/ZuK//JMjnqWZ8AEvUiku43fY6HPyTO6IcYEn
+                2Ume0/yeTnHkInfQ9Z2WqIdaacMthptuKHyXGvbNpZX4i7XZUlBOffb7dtJG0Qvyz+OSrQXbqaq/YR9x
+                Nm60NJowiBb8cY/cZHCN2/kxHfkJu+fqOzO4gBXLokkOOcsNPONou9fSXbn4fLyHX7iCDuxoXp/qE9Ti
+                Ir7pso7pgr7V+/pBt3T92BXtxH8H7Mmc''')
         with self.subTest('chemical-potential'):
             self.assertAlmostEqual64(args['η'], '''
-                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=''')
+                eNoBggF9/lU4VTgzOIo4jzhIOCc3Lzd5Ngw4xTf1N1o3YzV70II3ijZoNuIyZcpFybM1jzX0zgU28DVG
+                zyLKzMjGyXfICsjfx73Hl8dxMzgzgctEylLKFsqEyU3I2snuyJvHZMc7yLzHKctJy6rKqci4yGDIrcqK
+                yvLID8royADInMcIyJ7HI8jnx13HcMdax2fHctPiz7DJHzZENjI1V8jHx7AxHMpqx2PHpMgAyDk4SDgR
+                OH84gThhONA3ejY/OJw3Y84NyYY2C8x1NqU2BjfJyQXKYy9rN303ATZmNnY3TTaMNmI0e8mGyf7KKjYz
+                Nj80DDD+Mn3MhsqCM/bLsslNyRo3BzcuNUbLEMk7yDfI08nXyWjJu818zdTK+8hRyLPJt8jYx4/HEMiy
+                x5wzXTNbzP02CTdQNmPKFslCNc7MRcjSx5DJZ8hiOG04JDjJN4020jBZySAwQsvxyJDInMlnyCzIxMes
+                x3DHmMeAx1fHVMdfx0vHS8dsxzzHl8csxz3HSceixxfInGzEOA==''')

nutils/mesh.py

@@ -623,7 +623,7 @@ def unitsquare(nelems, etype):
     etype : :class:`str`
         Type of element used for meshing. Supported are:
 
-        * ``"square"``: structured mesh of squares.
+        * ``"square"`` or ``"rectilinear"``: structured mesh of squares.
 
         * ``"triangle"``: unstructured mesh of triangles.
 
@@ -652,7 +652,7 @@ def unitsquare(nelems, etype):
 
     space = 'X'
 
-    if etype == 'square':
+    if etype in ('square', 'rectilinear'):
         topo, geom = rectilinear([nelems, nelems], space=space)
         return topo, geom / nelems
 
@@ -711,11 +711,6 @@ def unitsquare(nelems, etype):
         raise Exception('invalid element type {!r}'.format(etype))
 
 
-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.
 
@@ -744,11 +739,42 @@ def unitcircle(nelems: int, variant: str) -> Tuple[Topology, function.Array]:
 
     if variant == 'rectilinear':
         topo, geom = unitsquare(nelems, 'square')
-        return topo, _square_to_circle(geom)
+        angle = (geom - 0.5) * (numpy.pi / 2)
+        return topo, numpy.sqrt(2) * numpy.sin(angle) * numpy.cos(angle)[[1, 0]]
 
     elif variant == 'multipatch':
         topo, geom = unitsquare(nelems, 'multipatch')
-        return topo, _square_to_circle(geom)
+
+        B, T, L, R, C = topo.basis('patch')
+        x, y = geom * 2 - 1
+
+        xlin = x / numpy.maximum(abs(y), 1/3) # -1 / 1
+        ylin = y / numpy.maximum(abs(x), 1/3) # -1 / 1
+        xcup = numpy.maximum(1.5 * abs(x) - .5, 0) # 1 \ 0 / 1
+        ycup = numpy.maximum(1.5 * abs(y) - .5, 0) # 1 \ 0 / 1
+
+        b = numpy.sqrt(1/3) # scales inner square
+        xx = (b + (1-b) * xcup)**2
+        yy = (b + (1-b) * ycup)**2
+
+        c = .5 * (numpy.sqrt(2) - 1) # scales outer radius
+        X = (R-L) * (xx + c * xcup**2 * (1 - ylin**2)) + (T+C+B) * xlin * yy
+        Y = (T-B) * (yy + c * ycup**2 * (1 - xlin**2)) + (L+C+R) * ylin * xx
+        W = 1 + c * (L+R) * xcup**2 * (1 + ylin**2) + c * (T+B) * ycup**2 * (1 + xlin**2)
+
+        # Rather than returning [X, Y] / W, we project the numerator and
+        # denominator onto a second order basis for efficient evaluation and
+        # correct boundary gradients at the patch interfaces.
+
+        basis = topo.basis('spline', 2)
+        # Minimize e = (f - B c)^2 = f^2 - 2 f B c + cT BT B c -> BT B c = BT f
+        BB, BX, BY, BW, XX, YY, WW = topo.integrate([basis[:,None] * basis[None,:],
+            basis * X, basis * Y, basis * W, X**2, Y**2, W**2], degree=4)
+        cx, cy, cw = [BB.solve(BF) for BF in (BX, BY, BW)]
+        ex, ey, ew = [FF + (BB @ cf - 2 * BF) @ cf for (cf, BF, FF) in ((cx, BX, XX), (cy, BY, YY), (cw, BW, WW))]
+        log.debug(f'NURBS projection errors: x={ex:.0e}, y={ey:.0e}, w={ew:.0e}')
+
+        return topo, (basis @ numpy.stack([cx, cy], 1)) / (basis @ cw)
 
     else:
         raise Exception('invalid variant {!r}'.format(variant))

tests/test_mesh.py

@@ -196,7 +196,7 @@ 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)
+        numpy.testing.assert_almost_equal(normalerr, 0, decimal=14)
 
 unitcircle(variant='rectilinear')
 unitcircle(variant='multipatch')