add multipatch topology

docs/conf.py

@@ -46,6 +46,7 @@ def __getattr__( self, attr ):
     'sphinx.ext.viewcode',
     'sphinx.ext.mathjax',
     'sphinx.ext.napoleon',
+    'sphinx.ext.intersphinx',
 ]
 
 # Add any paths that contain templates here, relative to this directory.
@@ -347,3 +348,9 @@ def __getattr__( self, attr ):
 
 autodoc_member_order = 'bysource'
 autodoc_default_flags = [ 'members' ]
+
+intersphinx_mapping = {
+    'python': ('https://docs.python.org/3', None),
+    'numpy': ('https://docs.scipy.org/doc/numpy/', None),
+    'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None)
+}

nutils/mesh.py

@@ -15,7 +15,7 @@
 """
 
 from . import topology, function, util, element, numpy, numeric, transform, log, _
-import os, warnings
+import os, warnings, itertools
 
 # MESH GENERATORS
 
@@ -95,6 +95,208 @@ def newrectilinear( nodes, periodic=None, bnames=[['left','right'],['bottom','to
     geom = function.concatenate( function.bifurcate(geom,geomi) )
   return domain, geom
 
+@log.title
+def multipatch( patches, nelems, patchverts=None, name='multipatch' ):
+  '''multipatch rectilinear mesh generator
+
+  Generator for a :class:`~nutils.topology.MultipatchTopology` and geometry.
+  The :class:`~nutils.topology.MultipatchTopology` consists of a set patches,
+  where each patch is a :class:`~nutils.topology.StructuredTopology` and all
+  patches have the same number of dimensions.
+
+  The ``patches`` argument, a :class:`numpy.ndarray`-like with shape
+  ``(npatches, 2*ndims)`` or ``(npatches,)+(2,)*ndims``, defines the
+  connectivity by labelling the patch vertices.  For example, three
+  one-dimensional patches can be connected at one edge by::
+
+      # connectivity:     3
+      #                   │
+      #                1──0──2
+
+      patches=[ [0,1], [0,2], [0,3] ]
+
+  Or two two-dimensional patches along an edge by::
+
+      # connectivity:  3──4──5
+      #                │  │  │
+      #                0──1──2
+
+      patches=[ [[0,3],[1,4]], [[1,4],[2,5]] ]
+
+  The geometry is specified by the ``patchverts`` argument: a
+  :class:`numpy.ndarray`-like with shape ``(nverts,ngeomdims)`` specifying for
+  each vertex a coordinate.  Note that the dimension of the geometry may be
+  higher than the dimension of the patches.  The created geometry is a
+  patch-wise linear interpolation of the vertex coordinates.  If the
+  ``patchverts`` argument is omitted the geometry describes a unit hypercube
+  per patch.
+
+  The ``nelems`` argument is either an :class:`int` defining the number of
+  elements per patch per dimension, or a :class:`dict` with edges (a pair of
+  vertex numbers) as keys and the number of elements (:class:`int`) as values,
+  with key ``None`` specifying the default number of elements.  Example::
+
+      # connectivity:  3─────4─────5
+      #                │ 4x3 │ 8x3 │
+      #                0─────1─────2
+
+      patches=[ [[0,3],[1,4]], [[1,4],[2,5]] ]
+      nelems={None: 4, (1,2): 8, (4,5): 8, (0,3): 3, (1,4): 3, (2,5): 3}
+
+  Since the patches are structured topologies, the number of elements per
+  patch per dimension should be unambiguous.  In above example specifying
+  ``nelems={None: 4, (1,2): 8}`` will raise an exception because the patch on
+  the right has 8 elements along edge ``(1,2)`` and 4 along ``(4,5)``.
+
+  Example
+  -------
+
+  An L-shaped domain can be generated by::
+
+      # connectivity:  2──5
+      #                │  |
+      #                1──4─────7     y
+      #                │  │     │     │
+      #                0──3─────6     └──x
+
+      domain, geom = mesh.multipatch(
+        patches=[ [0,1,3,4], [1,2,4,5], [3,4,6,7] ],
+        patchverts=[ [0,0], [0,1], [0,2], [1,0], [1,1], [1,2], [3,0], [3,1] ],
+        nelems={None: 4, (3,6): 8, (4,7): 8} )
+
+  The number of elements is chosen such that all elements in the domain have
+  the same size.
+
+  A topology and geometry describing the surface of a sphere can be generated
+  by creating a multipatch cube surface and inflating the cube to a sphere::
+
+      # connectivity:    3────7
+      #                 ╱│   ╱│
+      #                2────6 │     y
+      #                │ │  │ │     │
+      #                │ 1──│─5     │ z
+      #                │╱   │╱      │╱
+      #                0────4       *────x
+
+      topo, cube = multipatch(
+        patches=[
+          # The order of the vertices is chosen such that normals point outward.
+          [2,3,0,1],
+          [4,5,6,7],
+          [4,6,0,2],
+          [1,3,5,7],
+          [1,5,0,4],
+          [2,6,3,7],
+        ],
+        patchverts=tuple(itertools.product(*([[-1,1]]*3))),
+        nelems=10,
+      )
+      sphere = cube / function.sqrt((cube**2).sum(0))
+
+  Args
+  ----
+  patches:
+      A :class:`numpy.ndarray` with shape sequence of patches with each patch being a list of vertex indices.
+  patchverts:
+      A sequence of coordinates of the vertices.
+  nelems:
+      Either an :class:`int` specifying the number of elements per patch per
+      dimension, or a :class:`dict` with edges (a pair of vertex numbers) as
+      keys and the number of elements (:class:`int`) as values, with key
+      ``None`` specifying the default number of elements.
+
+  Returns
+  -------
+  :class:`nutils.topology.MultipatchTopology`:
+      The multipatch topology.
+  :class:`nutils.function.Array`:
+      The geometry defined by the ``patchverts`` or a unit hypercube per patch
+      if ``patchverts`` is not specified.
+  '''
+
+  patches = numpy.array( patches )
+  if patches.dtype != int:
+    raise ValueError( '`patches` should be an array of ints.' )
+  if len( patches.shape ) != 2 or patches.ndim > 2 and patches.shape[1:] != (2,)*(patches.ndim-1):
+    raise ValueError( '`patches` should be an array with shape (npatches,2,...,2) or (npatches,2*ndims).' )
+  patches = patches.reshape(patches.shape[0], -1)
+
+  # determine topological dimension of patches
+
+  ndims = 0
+  while 2**ndims < patches.shape[1]:
+    ndims += 1
+  if 2**ndims > patches.shape[1]:
+    raise ValueError( 'Only hyperrectangular patches are supported: ' \
+      'number of patch vertices should be a power of two.' )
+
+  # group all common patch edges (and/or boundaries?)
+
+  if isinstance( nelems, int ):
+    nelems = {None: nelems}
+  elif isinstance( nelems, dict ):
+    nelems = {(k and frozenset(k)): v for k, v in nelems.items()}
+  else:
+    raise ValueError( '`nelems` should be an `int` or `dict`')
+
+  # create patch topologies, geometries
+
+  if patchverts is not None:
+    patchverts = numpy.array( patchverts )
+    indices = set( patches.flat )
+    if tuple( sorted( indices ) ) != tuple( range( len( indices ) ) ):
+      raise ValueError( 'Patch vertices in `patches` should be numbered consecutively, starting at 0.' )
+    if len( patchverts ) != len( indices ):
+      raise ValueError( 'Number of `patchverts` does not equal number of vertices specified in `patches`.' )
+    if len( patchverts.shape ) != 2:
+      raise ValueError( 'Every patch vertex should be an array of dimension 1.' )
+
+  topos = []
+  coords = []
+  for i, patch in enumerate( patches ):
+    # find shape of patch and local patch coordinates
+    shaped_patch = patch.reshape( [2]*ndims )
+    shape = []
+    for dim in range( ndims ):
+      nelems_sides = []
+      sides = [(0,1)]*ndims
+      sides[dim] = slice(None),
+      for side in itertools.product(*sides):
+        sideverts = frozenset(shaped_patch[side])
+        if sideverts in nelems:
+          nelems_sides.append(nelems[sideverts])
+        else:
+          nelems_sides.append(nelems[None])
+      if len(set(nelems_sides)) != 1:
+        raise ValueError( 'duplicate number of elements specified for patch {} in dimension {}'.format( i, dim ) )
+      shape.append( nelems_sides[0] )
+    # create patch topology
+    topos.append( rectilinear( shape, name='{}{}'.format(name, i) )[0] )
+    # compute patch geometry
+    patchcoords = [ numpy.linspace(0, 1, n+1) for n in shape ]
+    patchcoords = numeric.meshgrid( *patchcoords ).reshape( ndims, -1 )
+    if patchverts is not None:
+      patchcoords = numpy.array([
+        sum(
+          patchverts[j]*util.product(c if s else 1-c for c, s in zip(coord, side))
+          for j, side in zip(patch, itertools.product(*[[0,1]]*ndims))
+        )
+        for coord in patchcoords.T
+      ]).T
+    coords.append( patchcoords )
+
+  # build patch boundary data
+
+  boundarydata = topology.MultipatchTopology.build_boundarydata( patch.reshape( (2,)*ndims ) for patch in patches )
+
+  # join patch topologies, geometries
+
+  topo = topology.MultipatchTopology( zip( topos, boundarydata ) )
+  funcsp = topo.splinefunc( degree=1, patchcontinuous=False )
+  geom = ( funcsp * numpy.concatenate( coords, axis=1 ) ).sum( -1 )
+
+  return topo, geom
+
 @log.title
 def gmsh( fname, name=None ):
   """Gmsh parser