From 9ff11cd2bbf3a7d3c1b8527378b3e064aca09bc7 Mon Sep 17 00:00:00 2001
From: Joost van Zwieten <joostvanzwieten@evalf.nl>
Date: Fri, 14 Apr 2017 14:13:08 +0200
Subject: [PATCH 1/5] fix TransformChain.lookup with empty transforms

---
 nutils/transform.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/nutils/transform.py b/nutils/transform.py
index 658280039..4fa6ec551 100644
--- a/nutils/transform.py
+++ b/nutils/transform.py
@@ -115,6 +115,8 @@ def promote( self, ndims ):
     return head, CanonicalTransformChain(tail)
 
   def lookup( self, transforms ):
+    if not transforms:
+      return
     for trans in transforms:
       ndims = trans.fromdims
       break

From 1d6fb53e11cf0c9469ac2219a0231af23ca65893 Mon Sep 17 00:00:00 2001
From: Joost van Zwieten <joostvanzwieten@evalf.nl>
Date: Sat, 15 Apr 2017 12:49:40 +0200
Subject: [PATCH 2/5] enable intersphinx, link to python, numpy, scipy

---
 docs/conf.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/docs/conf.py b/docs/conf.py
index 500b13fb8..fcf3c2c3d 100644
--- a/docs/conf.py
+++ b/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)
+}

From 6cd28edf92105253713d57ac21cd4accfb521e2a Mon Sep 17 00:00:00 2001
From: Joost van Zwieten <joostvanzwieten@evalf.nl>
Date: Fri, 14 Apr 2017 15:18:45 +0200
Subject: [PATCH 3/5] improve topology interface tests

---
 tests/topology.py | 40 +++++++++++++++++++++++-----------------
 1 file changed, 23 insertions(+), 17 deletions(-)

diff --git a/tests/topology.py b/tests/topology.py
index 4020e79c7..3640f3f67 100644
--- a/tests/topology.py
+++ b/tests/topology.py
@@ -53,17 +53,20 @@ def verify_boundaries( domain, geom ):
   numpy.testing.assert_array_almost_equal( lhs, rhs )
 
 def verify_interfaces( domain, geom, periodic ):
+  # If `periodic` is true, the domain should be a unit hypercube or this test
+  # might fail.  The function `f` defined below is C0 continuous on a periodic
+  # hypercube and Cinf continuous inside the hypercube.
   x1, x2, n1, n2 = domain.interfaces.elem_eval( [ geom, function.opposite(geom), geom.normal(), function.opposite(geom.normal()) ], 'gauss2', separate=False )
   if not periodic:
-    assert numpy.all( x1 == x2 )
-  assert numpy.all( n1 == -n2 )
+    numpy.testing.assert_array_almost_equal( x1, x2 )
+  numpy.testing.assert_array_almost_equal( n1, -n2 )
 
   # Test ∫_E f_,i = ∫_∂E f n_i ∀ E in `domain`.
   f = ((0.5 - geom)**2).sum(axis=0)
   elemindicator = domain.basis( 'discont', degree=0 ).vector( domain.ndims )
   lhs = domain.integrate( (elemindicator*f.grad(geom)[None]).sum(axis=1), ischeme='gauss2', geometry=geom )
   rhs = domain.interfaces.integrate( (-function.jump(elemindicator)*f*function.normal(geom)[None]).sum(axis=1), ischeme='gauss2', geometry=geom )
-  if not periodic:
+  if len( domain.boundary ):
     rhs += domain.boundary.integrate( (elemindicator*f*function.normal(geom)[None]).sum(axis=1), ischeme='gauss2', geometry=geom )
   numpy.testing.assert_array_almost_equal( lhs, rhs )
 
@@ -338,20 +341,23 @@ def test():
     numpy.testing.assert_array_almost_equal( located, target )
 
 
-@register( '3d', ndims=3 )
-@register( '2d_periodic', ndims=2, periodic=True )
-@register( '2d', ndims=2 )
-@register( '1d_periodic', ndims=1, periodic=True )
-#@register( '1d', ndims=1 ) # disabled, see issue #193
-def hierarchical( ndims, periodic=False ):
-
-  domain, geom = mesh.rectilinear( [numpy.linspace(0, 1, 5)]*ndims, periodic=range(ndims) if periodic else [] )
-  # Refine `domain`.
-  indicator = 1
-  for dim in range( ndims ):
-      indicator = function.min( indicator, 1 - geom[dim] )
-  for threshold in 0.5, 0.75:
-      domain = domain.refined_by( elem for elem, value in zip( domain, domain.elem_mean( [indicator], ischeme='gauss1', geometry=geom )[0] ) if value >= threshold )
+@register( '3d_l_rrr', pos=0, ndims=3 )
+@register( '3d_l_rpr', pos=0, ndims=3, periodic=[1] )
+@register( '2d_l_pp', pos=0, ndims=2, periodic=[0,1] )
+@register( '2d_l_pr', pos=0, ndims=2, periodic=[0] )
+@register( '2d_c_pr', pos=0.5, ndims=2, periodic=[0] )
+@register( '2d_r_pr', pos=1, ndims=2, periodic=[0] )
+@register( '2d_l_rr', pos=0, ndims=2 )
+@register( '1d_l_p', pos=0, ndims=1, periodic=[0] )
+@register( '1d_c_p', pos=0.5, ndims=1, periodic=[0] )
+#@register( '1d_l_r', pos=0, ndims=1 ) # disabled, see issue #193
+def hierarchical( ndims, pos, periodic=() ):
+
+  domain, geom = mesh.rectilinear( [numpy.linspace(0, 1, 7)]*ndims, periodic=periodic )
+  # Refine `domain` near `pos`.
+  distance = ((geom-pos)**2).sum(0)**0.5
+  for threshold in 0.3, 0.15:
+      domain = domain.refined_by( elem for elem, value in zip( domain, domain.elem_mean( [distance], ischeme='gauss1', geometry=geom )[0] ) if value <= threshold )
 
   if not periodic:
     @unittest

From 4cc18782c4709f10ec2613f6861b3042da471b82 Mon Sep 17 00:00:00 2001
From: Joost van Zwieten <joostvanzwieten@evalf.nl>
Date: Wed, 29 Jul 2015 19:49:39 +0200
Subject: [PATCH 4/5] add spline generator helper to StructuredTopology

---
 nutils/topology.py | 19 ++++++++++++-------
 1 file changed, 12 insertions(+), 7 deletions(-)

diff --git a/nutils/topology.py b/nutils/topology.py
index 704e050e7..7bdab2850 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -1079,8 +1079,8 @@ def interfaces( self ):
     assert len(itopos) == self.ndims
     return UnionTopology( itopos, names=[ 'dir{}'.format(idim) for idim in range(self.ndims) ] )
 
-  def basis_spline( self, degree, knotvalues=None, knotmultiplicities=None, periodic=None, removedofs=None ):
-    'spline basis'
+  def _basis_spline( self, degree, knotvalues=None, knotmultiplicities=None, periodic=None ):
+    'spline with structure information'
     
     if periodic is None:
       periodic = self.periodic
@@ -1096,11 +1096,6 @@ def basis_spline( self, degree, knotvalues=None, knotmultiplicities=None, period
     if knotmultiplicities is None:
       knotmultiplicities = [None]*self.ndims
 
-    if removedofs == None:
-      removedofs = [None] * self.ndims
-    else:
-      assert len(removedofs) == self.ndims
-
     vertex_structure = numpy.array( 0 )
     dofshape = []
     slices = []
@@ -1187,7 +1182,17 @@ def basis_spline( self, degree, knotvalues=None, knotmultiplicities=None, period
       dofs = vertex_structure[S].ravel()
       dofmap[trans] = dofs
       funcmap[trans] = std
+    return funcmap, dofmap, dofshape
+
+  def basis_spline( self, degree, knotvalues=None, knotmultiplicities=None, periodic=None, removedofs=None ):
+    'spline basis'
+
+    if removedofs == None:
+      removedofs = [None] * self.ndims
+    else:
+      assert len(removedofs) == self.ndims
 
+    funcmap, dofmap, dofshape = self._basis_spline( degree=degree, knotvalues=knotvalues, knotmultiplicities=knotmultiplicities, periodic=periodic )
     func = function.function( funcmap, dofmap, numpy.product(dofshape) )
     if not any( removedofs ):
       return func

From 5e3326370257f7e5bd35028bd0c90c59f4bcbefb Mon Sep 17 00:00:00 2001
From: Joost van Zwieten <joostvanzwieten@evalf.nl>
Date: Thu, 30 Jul 2015 18:02:49 +0200
Subject: [PATCH 5/5] add multipatch topology, mesh generator

---
 nutils/mesh.py     | 204 ++++++++++++++++++++++++++++++++++++++++++++-
 nutils/topology.py | 199 +++++++++++++++++++++++++++++++++++++++++++
 tests/topology.py  | 109 ++++++++++++++++++++++--
 3 files changed, 506 insertions(+), 6 deletions(-)

diff --git a/nutils/mesh.py b/nutils/mesh.py
index e4b3cb59b..025baf734 100644
--- a/nutils/mesh.py
+++ b/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
diff --git a/nutils/topology.py b/nutils/topology.py
index 7bdab2850..216459632 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -1917,6 +1917,205 @@ def __init__( self ):
   def basis( self, name, *args, **kwargs ):
     return function.asarray( [1] )
 
+class MultipatchTopology( Topology ):
+  'multipatch topology'
+
+  @staticmethod
+  def build_boundarydata( connectivity ):
+    'build boundary data based on connectivity'
+
+    boundarydata = []
+    for patch in connectivity:
+      ndims = len( patch.shape )
+      patchboundarydata = []
+      for dim, side in itertools.product( range( ndims ), [-1, 0] ):
+        # ignore vertices at opposite face
+        verts = numpy.array( patch )
+        opposite = tuple( {0:-1, -1:0}[side] if i == dim else slice(None) for i in range( ndims ) )
+        verts[opposite] = verts.max()+1
+        if len( set( verts.flat ) ) != 2**(ndims-1)+1:
+          raise NotImplementedError( 'Cannot compute canonical boundary if vertices are used more than once.' )
+        # reverse axes such that lowest vertex index is at first position
+        reverse = tuple( slice(None, None, -1) if i else slice(None) for i in numpy.unravel_index(verts.argmin(), verts.shape) )
+        verts = verts[reverse]
+        # transpose such that second lowest vertex connects to lowest vertex in first dimension, third in second dimension, et cetera
+        k = [ verts[tuple( 1 if i == j else 0 for j in range( ndims ) )] for i in range( ndims ) ]
+        transpose = tuple( sorted( range( ndims ), key=k.__getitem__ ) )
+        verts = verts.transpose( transpose )
+        # boundarid
+        boundaryid = tuple( verts[...,0].flat )
+        patchboundarydata.append( (boundaryid,dim,side,reverse,transpose) )
+      boundarydata.append( tuple( patchboundarydata ) )
+
+    # TODO: boundary sanity checks
+
+    return boundarydata
+
+  def __init__( self, patches ):
+    'constructor'
+
+    self.patches = tuple( patches )
+
+    self._patchinterfaces = {}
+    for topo, boundaries in self.patches:
+      for boundaryid, dim, side, reverse, transpose in boundaries:
+        self._patchinterfaces.setdefault( boundaryid, [] ).append(( topo, dim, side, reverse, transpose ))
+    self._patchinterfaces = {
+      boundaryid: tuple( data )
+      for boundaryid, data in self._patchinterfaces.items()
+      if len( data ) > 1
+    }
+
+    super().__init__( self.patches[0][0].ndims )
+
+  @cache.property
+  def elements( self ):
+    return tuple( itertools.chain.from_iterable( topo for topo, boundaries in self.patches ) )
+
+  def getitem( self, key ):
+    for i in range( len( self.patches ) ):
+      if key == 'patch{}'.format(i):
+        return self.patches[i][0]
+    else:
+      return UnionTopology( patch.getitem(key) for patch, boundaries in self.patches )
+
+  def basis_spline( self, degree, patchcontinuous=True ):
+    '''spline from vertices
+
+    Create a spline basis with degree ``degree`` per patch.  If
+    ``patchcontinuous``` is true the basis is $C^0$-continuous at patch
+    interfaces.
+    '''
+
+    funcmap = {}
+    dofmap = {}
+    dofcount = 0
+    commonboundarydofs = {}
+    for topo, boundaries in self.patches:
+      # build structured spline basis on patch `topo`
+      patchfuncmap, patchdofmap, patchdofcount = topo._basis_spline( degree )
+      funcmap.update( patchfuncmap )
+      # renumber dofs
+      dofmap.update( (trans,dofs+dofcount) for trans, dofs in patchdofmap.items() )
+      if patchcontinuous:
+        # reconstruct multidimensional dof structure
+        dofs = dofcount + numpy.arange( numpy.prod( patchdofcount ), dtype=int ).reshape( patchdofcount )
+        for boundaryid, dim, side, reverse, transpose in boundaries:
+          # get patch boundary dofs and reorder to canonical form
+          boundarydofs = dofs[reverse].transpose(transpose)[...,0].ravel()
+          # append boundary dofs to list (in increasing order, automatic by outer loop and dof increment)
+          commonboundarydofs.setdefault( boundaryid, [] ).append( boundarydofs )
+      dofcount += numpy.prod( patchdofcount )
+
+    if patchcontinuous:
+      # build merge mapping: merge common boundary dofs (from low to high)
+      pairs = itertools.chain(*( zip( *dofs ) for dofs in commonboundarydofs.values() if len( dofs ) > 1 ))
+      merge = {}
+      for dofs in sorted( pairs ):
+        dst = merge.get( dofs[0], dofs[0] )
+        for src in dofs[1:]:
+          merge[src] = dst
+      # build renumber mapping: renumber remaining dofs consecutively, starting at 0
+      remainder = set( merge.get( dof, dof ) for dof in range( dofcount ) )
+      renumber = dict( zip( sorted( remainder ), range( len( remainder ) ) ) )
+      # apply mappings
+      dofmap = dict(
+        (k, numpy.array( tuple( renumber[merge.get( dof, dof )] for dof in v.flat ), dtype=int ).reshape( v.shape ))
+        for k, v in dofmap.items() )
+      dofcount = len( remainder )
+
+    return function.function( funcmap, dofmap, dofcount )
+
+  def basis_discont( self, degree ):
+    'discontinuous shape functions'
+
+    funcmap = {}
+    dofmap = {}
+    dofcount = 0
+    for topo, boundaries in self.patches:
+      pbasis = topo.basis( 'discont', degree=degree )
+      funcmap.update( pbasis.func.stdmap )
+      dofmap.update( (trans,dofs+dofcount) for trans, dofs in pbasis.dofmap.dofmap.items() )
+      dofcount += len( pbasis )
+    return function.function( funcmap, dofmap, dofcount )
+
+  def basis_patch( self ):
+    'degree zero patchwise discontinuous basis'
+
+    fmap = {}
+    nmap = {}
+    for ipatch, ( topo, boundaries ) in enumerate( self.patches ):
+      fmap[ topo.root ] = util.product( element.PolyLine( element.PolyLine.bernstein_poly(0) ) for i in range( self.ndims ) )
+      nmap[ topo.root ] = numpy.array([ ipatch ])
+    return function.function( fmap, nmap, len( self.patches ) )
+
+  @cache.property
+  def boundary( self ):
+    'boundary'
+
+    subtopos = []
+    subnames = []
+    for i, (topo, boundaries) in enumerate( self.patches ):
+      names = dict( zip( itertools.product( range( self.ndims ), [0,-1] ), topo._bnames ) )
+      for boundaryid, dim, side, reverse, transpose in boundaries:
+        if boundaryid in self._patchinterfaces:
+          continue
+        subtopos.append( topo.boundary[names[dim,side]] )
+        subnames.append( 'patch{}-{}'.format( i, names[dim,side] ) )
+    if len( subtopos ) == 0:
+      return EmptyTopology( self.ndims-1 )
+    else:
+      return UnionTopology( subtopos, subnames )
+
+  @cache.property
+  def interfaces( self ):
+    '''interfaces
+
+    Return a topology with all element interfaces.  The patch interfaces are
+    accessible via the group ``'interpatch'`` and the interfaces *inside* a
+    patch via ``'intrapatch'``.
+    '''
+
+    intrapatchtopo = EmptyTopology( self.ndims-1 ) if not self.patches else \
+      UnionTopology( topo.interfaces for topo, boundaries in self.patches )
+
+    btopos = []
+    bconnectivity = []
+    for boundaryid, patchdata in self._patchinterfaces.items():
+      if len( patchdata ) > 2:
+        raise ValueError( 'Cannot create interfaces of multipatch topologies with more than two interface connections.' )
+      pairs = []
+      for topo, dim, side, reverse, transpose in patchdata:
+        names = dict( zip( itertools.product( range( self.ndims ), [0,-1] ), topo._bnames ) )
+        # get structured set of boundary elements
+        elems = topo.boundary[names[dim, side]].structure
+        # add singleton axis
+        elems = elems[tuple( _ if i == dim else slice( None ) for i in range( self.ndims ) )]
+        # apply canonical transformation
+        elems = elems[reverse].transpose(transpose)[..., 0]
+        shape = elems.shape
+        pairs.append( elems.flat )
+      # join element pairs
+      elems = [
+        element.Element( elem_a.reference, elem_a.transform, elem_b.transform, oriented=True )
+        for elem_a, elem_b in zip( *pairs )
+      ]
+      # create structured topology of joined element pairs
+      bpatch = numpy.array( boundaryid ).reshape( (2,)*(self.ndims-1) )
+      #btopos.append( StructuredTopology( numpy.array( elems ).reshape( shape ) ) )
+      btopos.append( UnstructuredTopology( self.ndims-1, elems ) )
+      bconnectivity.append( bpatch )
+    # create multipatch topology of interpatch boundaries
+    interpatchtopo = MultipatchTopology( zip( btopos, self.build_boundarydata( bconnectivity ) ) )
+
+    return UnionTopology( (intrapatchtopo, interpatchtopo), ('intrapatch', 'interpatch') )
+
+  @cache.property
+  def refined( self ):
+    'refine'
+
+    return MultipatchTopology( (topo.refined, boundaries) for topo, boundaries in self.patches )
+
 # UTILITY FUNCTIONS
 
 DimAxis = collections.namedtuple( 'DimAxis', ['i','j','isperiodic'] )
diff --git a/tests/topology.py b/tests/topology.py
index 3640f3f67..6ffda4905 100644
--- a/tests/topology.py
+++ b/tests/topology.py
@@ -1,7 +1,7 @@
 
 from nutils import *
 from . import register, unittest
-import numpy, copy, sys, pickle, subprocess, base64
+import numpy, copy, sys, pickle, subprocess, base64, itertools
 
 grid = numpy.linspace( 0., 1., 4 )
 
@@ -52,20 +52,24 @@ def verify_boundaries( domain, geom ):
   rhs = domain.boundary.integrate( f*function.normal(geom), ischeme='gauss2', geometry=geom )
   numpy.testing.assert_array_almost_equal( lhs, rhs )
 
-def verify_interfaces( domain, geom, periodic ):
+def verify_interfaces( domain, geom, periodic, interfaces=None, elemindicator=None ):
   # If `periodic` is true, the domain should be a unit hypercube or this test
   # might fail.  The function `f` defined below is C0 continuous on a periodic
   # hypercube and Cinf continuous inside the hypercube.
-  x1, x2, n1, n2 = domain.interfaces.elem_eval( [ geom, function.opposite(geom), geom.normal(), function.opposite(geom.normal()) ], 'gauss2', separate=False )
+  if interfaces is None:
+    interfaces = domain.interfaces
+  x1, x2, n1, n2 = interfaces.elem_eval( [ geom, function.opposite(geom), geom.normal(), function.opposite(geom.normal()) ], 'gauss2', separate=False )
   if not periodic:
     numpy.testing.assert_array_almost_equal( x1, x2 )
   numpy.testing.assert_array_almost_equal( n1, -n2 )
 
   # Test ∫_E f_,i = ∫_∂E f n_i ∀ E in `domain`.
   f = ((0.5 - geom)**2).sum(axis=0)
-  elemindicator = domain.basis( 'discont', degree=0 ).vector( domain.ndims )
+  if elemindicator is None:
+    elemindicator = domain.basis( 'discont', degree=0 )
+  elemindicator = elemindicator.vector( domain.ndims )
   lhs = domain.integrate( (elemindicator*f.grad(geom)[None]).sum(axis=1), ischeme='gauss2', geometry=geom )
-  rhs = domain.interfaces.integrate( (-function.jump(elemindicator)*f*function.normal(geom)[None]).sum(axis=1), ischeme='gauss2', geometry=geom )
+  rhs = interfaces.integrate( (-function.jump(elemindicator)*f*function.normal(geom)[None]).sum(axis=1), ischeme='gauss2', geometry=geom )
   if len( domain.boundary ):
     rhs += domain.boundary.integrate( (elemindicator*f*function.normal(geom)[None]).sum(axis=1), ischeme='gauss2', geometry=geom )
   numpy.testing.assert_array_almost_equal( lhs, rhs )
@@ -367,3 +371,98 @@ def boundaries():
   @unittest
   def interfaces():
     verify_interfaces( domain, geom, periodic )
+
+
+@register('3', npatches=(3,))
+@register('2x2', npatches=(2,2))
+@register('3x3', npatches=(3,3))
+@register('2x2x3', npatches=(2,2,3))
+def multipatch_hyperrect(npatches):
+
+  npatches = numpy.array(npatches)
+  indices = numpy.arange((npatches+1).prod()).reshape(npatches+1)
+
+  domain, geom = mesh.multipatch(
+    patches=[ indices[tuple(map(slice, i, numpy.array(i)+2))].ravel().tolist() for i in itertools.product(*map(range, npatches)) ],
+    patchverts=tuple( itertools.product( *map( range, npatches+1 ) ) ),
+    nelems=4,
+  )
+
+  @unittest
+  def spline_basis():
+    basis = domain.basis( 'spline', degree=2 )
+    coeffs = domain.project( 1, onto=basis, geometry=geom, ischeme='gauss4' )
+    numpy.testing.assert_array_almost_equal( coeffs, numpy.ones( coeffs.shape ) )
+
+  @unittest
+  def discont_basis():
+    basis = domain.basis( 'discont', degree=2 )
+    coeffs = domain.project( 1, onto=basis, geometry=geom, ischeme='gauss4' )
+    numpy.testing.assert_array_almost_equal( coeffs, numpy.ones( coeffs.shape ) )
+
+  @unittest
+  def boundaries():
+    verify_boundaries( domain, geom )
+
+  @unittest
+  def interfaces():
+    verify_interfaces( domain, geom, periodic=False )
+
+  @unittest
+  def interpatch_interfaces():
+    verify_interfaces( domain, geom, periodic=False, interfaces=domain.interfaces['interpatch'], elemindicator=domain.basis( 'patch' ) )
+
+
+@register
+def multipatch_L():
+
+  # 2---5
+  # |   |
+  # 1---4------7
+  # |   |      |
+  # 0---3------6
+
+  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} )
+
+  @unittest
+  def spline_basis():
+    basis = domain.basis( 'spline', degree=2 )
+    coeffs = domain.project( 1, onto=basis, geometry=geom, ischeme='gauss4' )
+    numpy.testing.assert_array_almost_equal( coeffs, numpy.ones( coeffs.shape ) )
+
+  @unittest
+  def discont_basis():
+    basis = domain.basis( 'discont', degree=2 )
+    coeffs = domain.project( 1, onto=basis, geometry=geom, ischeme='gauss4' )
+    numpy.testing.assert_array_almost_equal( coeffs, numpy.ones( coeffs.shape ) )
+
+  @unittest
+  def patch_basis():
+    patch_index = domain.basis( 'patch' ).dot([0, 1, 2])
+    for ipatch in range(3):
+        vals = domain['patch{}'.format(ipatch)].elem_eval( patch_index, ischeme='gauss1' )
+        numpy.testing.assert_array_almost_equal( vals, ipatch )
+
+
+#@register  # not yet implemented
+def multipatch_corner_connectivity():
+
+  #     4---6
+  #     |   |
+  # 1---3---5
+  # |   |
+  # 0---2
+
+  domain, geom = mesh.multipatch(
+    patches=[ [0,1,2,3], [3,4,5,6] ],
+    patchverts=[ [0,0], [0,1], [1,0], [1,1], [1,2], [2,1], [2,2] ],
+    nelems=1 )
+
+  @unittest
+  def test():
+    basis = domain.basis('spline', degree=1)
+    # 4 basis functions per patch, minus one for vertex 3.
+    assert len(basis) == 7