diff --git a/examples/drivencavity-compatible.py b/examples/drivencavity-compatible.py
deleted file mode 100644
index 1a2e9c940..000000000
--- a/examples/drivencavity-compatible.py
+++ /dev/null
@@ -1,157 +0,0 @@
-#! /usr/bin/env python3
-#
-# In this script we solve the lid driven cavity problem for stationary Stokes
-# and Navier-Stokes flow. That is, a unit square domain, with no-slip left,
-# bottom and right boundaries and a top boundary that is moving at unit
-# velocity in positive x-direction.
-#
-# The script is identical to :ref:`examples/drivencavity.py` except that it
-# uses the Raviart-Thomas discretization providing compatible velocity and
-# pressure spaces resulting in a pointwise divergence-free velocity field.
-
-from nutils import mesh, function, solver, export, cli, testing
-from nutils.expression_v2 import Namespace
-import numpy
-import treelog
-
-
-def main(nelems: int, degree: int, reynolds: float):
-    '''
-    Driven cavity benchmark problem using compatible spaces.
-
-    .. arguments::
-
-       nelems [12]
-         Number of elements along edge.
-       degree [2]
-         Polynomial degree for velocity; the pressure space is one degree less.
-       reynolds [1000]
-         Reynolds number, taking the domain size as characteristic length.
-    '''
-
-    verts = numpy.linspace(0, 1, nelems+1)
-    domain, geom = mesh.rectilinear([verts, verts])
-
-    ns = Namespace()
-    ns.δ = function.eye(domain.ndims)
-    ns.Σ = function.ones([domain.ndims])
-    ns.x = geom
-    ns.Re = reynolds
-
-    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-    ns.add_field(('u', 'v'), function.vectorize([
-        domain.basis('spline', degree=(degree, degree-1), removedofs=((0, -1), None)),
-        domain.basis('spline', degree=(degree-1, degree), removedofs=(None, (0, -1)))]))
-    ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1))
-    ns.add_field(('λ', 'μ'))
-
-    ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
-    ns.uwall = numpy.stack([domain.boundary.indicator('top'), 0])
-    ns.N = 5 * degree * nelems  # nitsche constant based on element size = 1/nelems
-    ns.nitsche_i = '(N v_i - (∇_j(v_i) + ∇_i(v_j)) n_j) / Re'
-
-    res = domain.integral('∇_j(v_i) σ_ij dV' @ ns, degree=2*degree)
-    res += domain.boundary.integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dS' @ ns, degree=2*degree)
-    res += domain.integral('q (∇_k(u_k) + λ) dV' @ ns, degree=2*degree)
-    res += domain.integral('μ p dV' @ ns, degree=2*degree)
-
-    with treelog.context('stokes'):
-        args0 = solver.solve_linear('u:v,p:q,λ:μ', res)
-        postprocess(domain, ns, **args0)
-
-    res += domain.integral('v_i ∇_j(u_i) u_j dV' @ ns, degree=3*degree)
-    with treelog.context('navierstokes'):
-        args1 = solver.newton('u:v,p:q,λ:μ', res, arguments=args0).solve(tol=1e-10)
-        postprocess(domain, ns, **args1)
-
-    return args0, args1
-
-# Postprocessing in this script is separated so that it can be reused for the
-# results of Stokes and Navier-Stokes, and because of the extra steps required
-# for establishing streamlines.
-
-
-def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
-
-    div = domain.integral('∇_k(u_k)^2 dV' @ ns, degree=1).eval(**arguments)**.5
-    treelog.info('velocity divergence: {:.2e}'.format(div))  # confirm that velocity is pointwise divergence-free
-
-    ns = ns.copy_() # copy namespace so that we don't modify the calling argument
-    ns.add_field('stream', domain.basis('std', degree=2)[1:]) # stream function (sans first dof to obtain non-singular system)
-    ns.ε = function.levicivita(2)
-    sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(stream))^2 dV' @ ns, degree=4)
-    arguments = solver.optimize('stream,', sqr, arguments=arguments) # compute streamlines
-
-    bezier = domain.sample('bezier', 9)
-    x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'stream'] @ ns, **arguments)
-    with export.mplfigure('flow.png') as fig:  # plot velocity as field, pressure as contours, streamlines as dashed
-        ax = fig.add_axes([.1, .1, .8, .8], yticks=[], aspect='equal')
-        import matplotlib.collections
-        ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2))
-        ax.tricontour(x[:, 0], x[:, 1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9)
-        caxu = fig.add_axes([.1, .1, .03, .8], title='velocity')
-        imu = ax.tripcolor(x[:, 0], x[:, 1], bezier.tri, u, shading='gouraud', cmap='jet')
-        fig.colorbar(imu, cax=caxu)
-        caxu.yaxis.set_ticks_position('left')
-        caxp = fig.add_axes([.87, .1, .03, .8], title='pressure')
-        imp = ax.tricontour(x[:, 0], x[:, 1], bezier.tri, p, 16, cmap='gray', linestyles='solid')
-        fig.colorbar(imp, cax=caxp)
-
-# If the script is executed (as opposed to imported), :func:`nutils.cli.run`
-# calls the main function with arguments provided from the command line. To
-# keep with the default arguments simply run :sh:`python3
-# drivencavity-compatible.py`.
-
-
-if __name__ == '__main__':
-    cli.run(main)
-
-# Once a simulation is developed and tested, it is good practice to save a few
-# strategic return values for regression testing. The :mod:`nutils.testing`
-# module, which builds on the standard :mod:`unittest` framework, facilitates
-# this by providing :func:`nutils.testing.TestCase.assertAlmostEqual64` for the
-# embedding of desired results as compressed base64 data.
-
-
-class test(testing.TestCase):
-
-    def test_p2(self):
-        args0, args1 = main(nelems=3, reynolds=100, degree=2)
-        with self.subTest('stokes-velocity'):
-            self.assertAlmostEqual64(args0['u'], '''
-                eNrzu9Bt8OuUndkD/eTTSqezzP2g/E3698/ZmZlf2GjSaHJS3/90/Wm/C4qGh04CAErOF+g=''')
-        with self.subTest('stokes-pressure'):
-            self.assertAlmostEqual64(args0['p'], '''
-                eNoBIADf/0fSMC4P0ZLKjCk4JC7Pwy901sjb0jA90Lkt0NHxLm41+KgP+Q==''')
-        with self.subTest('stokes-multiplier'):
-            self.assertAlmostEqual(args0['λ'], 0)
-        with self.subTest('navier-stokes-velocity'):
-            self.assertAlmostEqual64(args1['u'], '''
-                eNoBMADP/2XOWjJSy5k1jS+yyzvLODfgL1rO0MrINpsxHM2ZNSrPqDTANCPVQsxCzeAvcc04yT3AF9c=''')
-        with self.subTest('navier-stokes-pressure'):
-            self.assertAlmostEqual64(args1['p'], '''
-                eNpbqpasrWa46TSTfoT+krO/de/pntA3Pnf2zFNtn2u2hglmAOKVDlE=''')
-        with self.subTest('navier-stokes-multiplier'):
-            self.assertAlmostEqual(args1['λ'], 0)
-
-    def test_p3(self):
-        args0, args1 = main(nelems=3, reynolds=100, degree=3)
-        with self.subTest('stokes-velocity'):
-            self.assertAlmostEqual64(args0['u'], '''
-                eNqbo3f7/AeDb2dmGGtd7r+odk7icoQJgjUHLpty8b/BvrMzjF/rll9qMD5kwAAFopc6dRvO2J2fo8d4
-                3sko4wwAjL4lyw==''', atol=1e-14)
-        with self.subTest('stokes-pressure'):
-            self.assertAlmostEqual64(args0['p'], '''
-                eNpbrN+us+Z8qWHMyfcXrly013l1ZocRAxwI6uvoHbwsZuxxNvZC5eUQg+5zS8wAElAT9w==''')
-        with self.subTest('stokes-multiplier'):
-            self.assertAlmostEqual(args0['λ'], 0)
-        with self.subTest('navier-stokes-velocity'):
-            self.assertAlmostEqual64(args1['u'], '''
-                eNoBUACv/14yGcxyNPbJYTahLj/LSDE7yy43SM9WMsXJoDR+N3Iw8s1hM5zJODeizcE0X8phNrQwUDOO
-                NbMzJi+ty4s1oDFqzxIzysjWzXIwFM3tNMjIpDMlJw==''')
-        with self.subTest('navier-stokes-pressure'):
-            self.assertAlmostEqual64(args1['p'], '''
-                eNoBMgDN/yoxvDHizaEy88kDLgoviCt4znIzMy7jL7nTRMzqzhAwBTA/LE0w6czY2GQwoNEYMHE3NDUW
-                DQ==''')
-        with self.subTest('navier-stokes-multiplier'):
-            self.assertAlmostEqual(args1['λ'], 0)
diff --git a/examples/drivencavity.py b/examples/drivencavity.py
index fc2731b29..ea7ac718a 100644
--- a/examples/drivencavity.py
+++ b/examples/drivencavity.py
@@ -1,23 +1,73 @@
 #! /usr/bin/env python3
 #
 # In this script we solve the lid driven cavity problem for stationary Stokes
-# and Navier-Stokes flow. That is, a unit square domain, with no-slip left,
-# bottom and right boundaries and a top boundary that is moving at unit
-# velocity in positive x-direction.
+# and Navier-Stokes flow. This benchmark problem consists of a square domain
+# with fixed left, bottom and right boundaries and a top boundary that is
+# moving at unit velocity in positive x-direction. Reference results can be
+# found for instance at https://www.acenumerics.com/the-benchmarks.html.
+#
+# The general conservation laws are:
+#
+#     Dρ/Dt = 0 (mass)
+#     ρ Du_i/Dt = ∇_j(σ_ij) (momentum)
+#
+# where we used the material derivative D·/Dt := ∂·/∂t + ∇_j(· u_j). The stress
+# tensor is σ_ij := μ (∇_j(u_i) + ∇_i(u_j) - 2 ∇_k(u_k) δ_ij / δ_nn) - p δ_ij,
+# with pressure p and dynamic viscosity μ, and ρ is the density.
+#
+# Assuming steady, incompressible flow, we take density to be constant. Further
+# introducing a reference length L and reference velocity U, we make the
+# equations dimensionless by taking spatial coordinates relative to L, velocity
+# relative to U, and pressure relative to ρ U^2. This reduces the conservation
+# laws to:
+#
+#     ∇_k(u_k) = 0 (mass)
+#     Du_i/Dt = ∇_j(σ_ij) (momentum)
+#
+# where the material derivative simplifies to D·/Dt := ∇_j(·) u_j, and the
+# stress tensor becomes σ_ij := (∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij, with
+# Reynolds number Re = ρ U L / μ.
+#
+# The weak form is obtained by multiplication of a test function and partial
+# integration of the right hand side of the momentum balance.
+#
+#     ∀ q: ∫_Ω q ∇_k(u_k) = 0 (mass)
+#     ∀ v: ∫_Ω (v_i Du_i/Dt + ∇_j(v_i) σ_ij) = ∫_Γ v_i σ_ij n_j (momentum)
+#
+# A remaining issue with this system is that its solution is not unique due to
+# its strictly kinematic boundary conditions: since ∫_Ω ∇_k(u_k) = 0 for any u
+# that satisfies the non-penetrating boundary conditions, any pressure space
+# that contains unity results in linear dependence. Furthermore, the strong
+# form shows that constant pressure shifts do not affect the momentum balance.
+# Both issues are solved by arbitrarily removing one of the basis functions.
+#
+# The normal velocity components at the wall are strongly constrained. The
+# tangential components are either strongly or weakly constrained depending on
+# the `strongbc` parameter. Since the tangential components at the top are
+# incompatible with the normal components at the left and right boundary, the
+# constraints are constructed in two steps to avoid Gibbs type oscillations,
+# and to make sure that the non-penetrating condition takes precedence.
+#
+# Depending on the `compatible` parameter, the script uses either a Taylor-Hood
+# (False) or a Raviart-Thomas (True) discretization. In case of TH, the system
+# is consistently modified by adding ∫_Ω .5 u_i v_i ∇_j(u_j) to yield the skew-
+# symmetric advective term ∫_Ω .5 (v_i ∇_j(u_i) - u_i ∇_j(v_i)) u_j. In case of
+# RT, the discretization guarantees a pointwise divergence-free velocity field,
+# and skew symmetry is implied.
 
 from nutils import mesh, function, solver, export, cli, testing
 from nutils.expression_v2 import Namespace
+import treelog as log
 import numpy
-import treelog
 
 
-def main(nelems: int, etype: str, degree: int, reynolds: float):
+def main(nelems: int, etype: str, degree: int, reynolds: float, compatible: bool, strongbc: bool):
     '''
     Driven cavity benchmark problem.
 
     .. arguments::
 
-       nelems [12]
+       nelems [32]
          Number of elements along edge.
        etype [square]
          Element type (square/triangle/mixed).
@@ -25,36 +75,69 @@ def main(nelems: int, etype: str, degree: int, reynolds: float):
          Polynomial degree for velocity; the pressure space is one degree less.
        reynolds [1000]
          Reynolds number, taking the domain size as characteristic length.
+       strongbc [no]
+         Use strong boundary constraints
+       compatible [no]
+         Use compatible spaces and weakly imposed boundary conditions; requires
+         etype=square and strongbc=no
     '''
 
+    if compatible and (strongbc or etype != 'square'):
+        raise Exception(f'compatible mode requires square elements and weak boundary conditions')
+
     domain, geom = mesh.unitsquare(nelems, etype)
 
     ns = Namespace()
     ns.δ = function.eye(domain.ndims)
+    ns.ε = function.levicivita(2)
     ns.Σ = function.ones([domain.ndims])
     ns.Re = reynolds
+    ns.uwall = numpy.stack([domain.boundary.indicator('top'), 0])
     ns.x = geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-    ns.add_field(('u', 'v'), domain.basis('std', degree=degree), shape=(domain.ndims,))
-    ns.add_field(('p', 'q'), domain.basis('std', degree=degree-1))
+    if not compatible:
+        ns.add_field(('u', 'v'), domain.basis('std', degree=degree), shape=(domain.ndims,))
+        ns.add_field(('p', 'q'), domain.basis('std', degree=degree-1)[1:])
+        ns.add_field('ψ', domain.basis('std', degree=2)[1:])
+    else:
+        ns.add_field(('u', 'v'), function.vectorize([
+            domain.basis('spline', degree=(degree, degree-1)),
+            domain.basis('spline', degree=(degree-1, degree))]))
+        ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1)[1:])
+        ns.add_field('ψ', domain.basis('spline', degree=degree)[1:])
     ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
 
-    sqr = domain.boundary.integral('u_k u_k dS' @ ns, degree=degree*2)
-    sqr += domain.locate(ns.x, [(1,1)], weights=numpy.array([1.]), tol=1e-10).integral(ns.p**2)
-    cons = solver.optimize('u,p', sqr, droptol=1e-15)
+    # weak formulation for Stokes flow, over-integrating for improved
+    # efficiency when we turn this to Navier-Stokes later on
+    res = domain.integral('∇_j(v_i) σ_ij dV' @ ns, degree=degree*3)
+    res += domain.integral('q ∇_k(u_k) dV' @ ns, degree=degree*3)
 
-    # project lid velocity separately to avoid Gibbs wiggles around boundary incompatibilities
-    sqr = domain.boundary['top'].integral('(u_0 - 1)^2 dS' @ ns, degree=degree*2)
-    lidcons = solver.optimize('u,', sqr, droptol=1e-15)
-    cons['u'] = numpy.choose(numpy.isnan(lidcons['u']), [lidcons['u'], cons['u']])
+    # strong enforcement of non-penetrating boundary conditions
+    sqr = domain.boundary.integral('(u_k n_k)^2 dS' @ ns, degree=degree*2)
+    cons = solver.optimize('u,', sqr, droptol=1e-15)
 
-    res = domain.integral('(∇_j(v_i) σ_ij + q ∇_k(u_k)) dV' @ ns, degree=degree*2)
-    with treelog.context('stokes'):
+    if strongbc:
+        # strong enforcement of tangential boundary conditions
+        sqr = domain.boundary.integral('(ε_ij n_i (u_j - uwall_j))^2 dS' @ ns, degree=degree*2)
+        tcons = solver.optimize('u,', sqr, droptol=1e-15)
+        cons['u'] = numpy.choose(numpy.isnan(cons['u']), [cons['u'], tcons['u']])
+    else:
+        # weak enforcement of tangential boundary conditions via Nitsche's method
+        ns.N = 5 * degree * nelems # Nitsche constant based on element size = 1/nelems
+        ns.nitsche_i = '(N v_i - (∇_j(v_i) + ∇_i(v_j)) n_j) / Re'
+        res += domain.boundary.integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dS' @ ns, degree=2*degree)
+
+    with log.context('stokes'):
         args0 = solver.solve_linear('u:v,p:q', res, constrain=cons)
         postprocess(domain, ns, **args0)
 
-    res += domain.integral('.5 (v_i ∇_j(u_i) - u_i ∇_j(v_i)) u_j dV' @ ns, degree=degree*3)
-    with treelog.context('navierstokes'):
+    # change to Navier-Stokes by adding convection
+    res += domain.integral('v_i ∇_j(u_i) u_j dV' @ ns, degree=degree*3)
+    if not compatible:
+        # add consistent term for skew-symmetry
+        res += domain.integral('.5 u_i v_i ∇_j(u_j) dV' @ ns, degree=degree*3)
+
+    with log.context('navier-stokes'):
         args1 = solver.newton('u:v,p:q', res, arguments=args0, constrain=cons).solve(tol=1e-10)
         postprocess(domain, ns, **args1)
 
@@ -64,35 +147,43 @@ def main(nelems: int, etype: str, degree: int, reynolds: float):
 # results of Stokes and Navier-Stokes, and because of the extra steps required
 # for establishing streamlines.
 
+def postprocess(domain, ns, **arguments):
 
-def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
-
-    ns = ns.copy_() # copy namespace so that we don't modify the calling argument
-    ns.add_field('stream', domain.basis('std', degree=2)[1:]) # stream function (sans first dof to obtain non-singular system)
-    ns.ε = function.levicivita(2)
-    sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(stream))^2 dV' @ ns, degree=4)
-    arguments = solver.optimize('stream,', sqr, arguments=arguments) # compute streamlines
+    # reconstruct velocity streamlines
+    sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(ψ))^2 dV' @ ns, degree=4)
+    arguments = solver.optimize('ψ,', sqr, arguments=arguments)
 
     bezier = domain.sample('bezier', 9)
-    x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'stream'] @ ns, **arguments)
-    with export.mplfigure('flow.png') as fig:  # plot velocity as field, pressure as contours, streamlines as dashed
-        ax = fig.add_axes([.1, .1, .8, .8], yticks=[], aspect='equal')
-        import matplotlib.collections
-        ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2))
-        ax.tricontour(x[:, 0], x[:, 1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9)
-        caxu = fig.add_axes([.1, .1, .03, .8], title='velocity')
-        imu = ax.tripcolor(x[:, 0], x[:, 1], bezier.tri, u, shading='gouraud', cmap='jet')
-        fig.colorbar(imu, cax=caxu)
-        caxu.yaxis.set_ticks_position('left')
-        caxp = fig.add_axes([.87, .1, .03, .8], title='pressure')
-        imp = ax.tricontour(x[:, 0], x[:, 1], bezier.tri, p, 16, cmap='gray', linestyles='solid')
-        fig.colorbar(imp, cax=caxp)
+    x, u, p, ψ = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'ψ'] @ ns, **arguments)
+
+    with export.mplfigure('flow.png', dpi=150) as fig: # plot velocity as field, pressure as contours, streamlines as dashed
+        ax = fig.add_subplot(111, aspect='equal')
+        im = export.triplot(ax, x, u, tri=bezier.tri, hull=bezier.hull, cmap='hot_r', clim=(0,1))
+        fig.colorbar(im, label='velocity')
+        ax.tricontour(x[:, 0], x[:, 1], bezier.tri, ψ, levels=numpy.percentile(ψ, numpy.arange(2,100,3)), colors='k', linestyles='solid', linewidths=.5, zorder=9)
+
+    x = numpy.linspace(0, 1, 1001)
+    v = domain.locate(ns.x, numpy.stack([x, numpy.repeat(.5, len(x))], 1), tol=1e-10).eval(ns.u[1], **arguments)
+    imin = numpy.argmin(v)
+    imax = numpy.argmax(v)
+    with export.mplfigure('cross-hor.png', dpi=150) as fig:
+        ax = fig.add_subplot(111, xlim=(0,1), title='horizontal cross section at y=0.5', xlabel='x-coordinate', ylabel='vertical velocity')
+        ax.plot(x, v)
+        ax.annotate(f'x={x[imax]:.3f}\nv={v[imax]:.3f}', xy=(x[imax], v[imax]), xytext=(x[imax], 0), arrowprops=dict(arrowstyle="->"), ha='center', va='center')
+        ax.annotate(f'x={x[imin]:.3f}\nv={v[imin]:.3f}', xy=(x[imin], v[imin]), xytext=(x[imin], 0), arrowprops=dict(arrowstyle="->"), ha='center', va='center')
+
+    y = numpy.linspace(0, 1, 1001)
+    u = domain.locate(ns.x, numpy.stack([numpy.repeat(.5, len(y)), y], 1), tol=1e-10).eval(ns.u[0], **arguments)
+    imin = numpy.argmin(u)
+    with export.mplfigure('cross-ver.png', dpi=150) as fig:
+        ax = fig.add_subplot(111, ylim=(0,1), title='vertical cross section at x=0.5', ylabel='y-coordinate', xlabel='horizontal velocity')
+        ax.plot(u, y)
+        ax.annotate(f'y={y[imin]:.3f}\nu={u[imin]:.3f}', xy=(u[imin], y[imin]), xytext=(0, y[imin]), arrowprops=dict(arrowstyle="->"), ha='left', va='center')
 
 # If the script is executed (as opposed to imported), :func:`nutils.cli.run`
 # calls the main function with arguments provided from the command line. To
 # keep with the default arguments simply run :sh:`python3 drivencavity.py`.
 
-
 if __name__ == '__main__':
     cli.run(main)
 
@@ -102,50 +193,82 @@ def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
 # this by providing :func:`nutils.testing.TestCase.assertAlmostEqual64` for the
 # embedding of desired results as compressed base64 data.
 
-
 class test(testing.TestCase):
 
-    def test_square(self):
-        args0, args1 = main(nelems=3, etype='square', reynolds=100, degree=3)
+    def test_baseline(self):
+        args0, args1 = main(nelems=3, etype='square', degree=2, reynolds=100, compatible=False, strongbc=False)
         with self.subTest('stokes-velocity'):
             self.assertAlmostEqual64(args0['u'], '''
-                eNpjYCAMgswhdJ1O+uWtl7/rp13hMyq4rmx8/cI/44qLW0wMLliZMhrxGMHUPT//WmfruW8GWudSjJ6d
-                FTeWP/vf+NH5TuNVhurGPqZvL8LUbTi3X+P1WV2D2rPthjZnw4yenflpdODSFaO9RpuNZple14Wp+362
-                VzP87Ce9b2f0DHaduW+w7EyEIYPeH4MW4z6Dh6apSOqKr4Wf5bv47cyl87vOKJ5fdmbFOQY9lvMtxkXn
-                H5rOvoSw1/H667OXz9eerTxnc3bV2Wdn2M8euKRzdq+R79lZppqXEP4Qvbz1HNd5rXNzzj47+/KM/FnG
-                M4/Ol59ZZXjzjI+psB4iXGbqbL3MeSHtyqezBdfvnrl+gelMxUWf0wYXjp1iNPpyFqaOmHAGAFFhkvE=''')
+                eNo9jSsLwlAcxQ82gwOjoFZ1ONA9NK6vC2Ky2cU2QeziFxDBJIKfZLv/O92DYRqo6waDGBR3GXjKgfPg
+                Bwh9A6AYA80TMOwCU0OkcwVonGMl8Uaa5C90i9+MipP2tn1gEgFXvgpDGqgRK+smPTWbaqqcdWMfeBPw
+                YcDGBdLMJR3Y/X+zdkhHHrEHM6lENt25+OU8OUi8PUn+klm87lacqiN4uQrZ4tUCLh3g4AFrV6Q/uctG
+                gQ==''')
         with self.subTest('stokes-pressure'):
             self.assertAlmostEqual64(args0['p'], '''
-                eNp7dOb9mRdnHp2pPbPw9MzTN848OXPpzJ4z1Wcizqw/c//Ma6DM9TMHzsw/s+PMFxTIdfbvGfazwmf1
-                zkaftTgrdJblrORZ47NTz94463qW/ezPM4xAcsLZjZcZGAAX7kL6''')
+                eNoBHgDh/+vVsNEXy6jTbiq1z7Av9C0mLJkw1NDTLEEtEC/xNAwED0s=''')
         with self.subTest('navier-stokes-velocity'):
             self.assertAlmostEqual64(args1['u'], '''
-                eNpjYCAMgswh9ErtA5c9ruTp/7xiZbhX28jo8MVbRn1XLIxVLhcbBxuIGsDUHbhgoxNx3tnA9vwcQ4fz
-                PkbC59mNP2rONOIxnGiUaOx9GaYu5PxOjfJzB/Xdz5kbWp9jNYo9t81ont4so1OGXUbNJtX6MHVd515p
-                XT4rqT//7EWDprPzDR+eTTY6pBdltNhoq9ELE+3zMHWe58rVl53lv2R1Vvbq+zOTdCLPmhpONkgwlDMO
-                NawyvWAIU/f7nMRN03PyF3rOSp6XOqt3bv+ZA+cirnSclzf2vxhvGgo3T/pC5xXW80Ln751ddzb1rOpZ
-                wzOXTp89l3iu1vic0WrTImOYugCd8uvrriy/aHf1w9nkizPPvDl/7rTe+cRTBmf3nVQx0T4DU0dMOAMA
-                p1CDeg==''')
+                eNpjYAABz4sMDP81GRh69BgYTgBxkQlIdCVQ5Nb5R/oSV+KNW8/+NXlidMOE4WS+mZM5A8M1JQaGuLMt
+                ejPOpBgtOs1vonRe0IT59ErjfWZAuasMDB1n9M9vOG2kL3B6l1H0uUXGRQZXjVYB9a2/xMCQfW7FueVn
+                nM/5nf5xevqZxDOq5w4bCwLlOoD6XDV/n1t//s5ZvjPzTjmdDjx55+Slky/MGaDgHFB3vz4DgynQfS9O
+                A3WcBIkCAB7aSkk=''')
         with self.subTest('navier-stokes-pressure'):
             self.assertAlmostEqual64(args1['p'], '''
-                eNoNiT0OQEAYBdd5xD1cSOIMiIJGVJuIiliS9VNYKolkNG6j9BUvmZlnmJnoSAnx6RmFNSUVjgErRVOQ
-                iFtWKTUZMQ2n/BuGnJaFi52bl48D73Fis+wjCpT6AWgsRHE=''')
+                eNpz01W9oHVmuU7SJYtzgherdcr0n59dfiZT11yP97yCGQDN0Azu''')
 
     def test_mixed(self):
-        args0, args1 = main(nelems=3, etype='mixed', reynolds=100, degree=2)
+        args0, args1 = main(nelems=3, etype='mixed', degree=2, reynolds=100, compatible=False, strongbc=False)
+        with self.subTest('stokes-velocity'):
+            self.assertAlmostEqual64(args0['u'], '''
+                eNpjYAABHx0Ghrbz+lcYGJZpR2hrnDtm/EObgWHvWSGD1WeuGTIwmJ9jYLAwnH32p8nKMzFmFqfVTJTP
+                aBszMOw0AenebmJh9tuMgWHGWX9jQ3MGhr9nLA35zxRcWHOm4mzOBQaG55cZGCTPGV6fdUrhwtEzvhe2
+                n+Y8k3RG+Mwio99nuoHqg4G48WzCmTignYUXDfXNzoedATuL4bMeA0Op9qczWqfXnTl2ioHhINAdHufv
+                ntx18qoZSH7FSRAJAB13Sc0=''')
+        with self.subTest('stokes-pressure'):
+            self.assertAlmostEqual64(args0['p'], '''
+                eNp7pKl+nf1KznmxS62ns/W+az/TNTL4ondU1/46t6GKKQDiJg1H''')
+        with self.subTest('navier-stokes-velocity'):
+            self.assertAlmostEqual64(args1['u'], '''
+                eNpjYACBVVcYGH5fZAKSChcLLv05p2jsp8XA4Hhe4YrV2QmGDAxHzzMwuBpknHcz4T3nYGp7JslY+ewy
+                YwaGamOQ7jyjzaZJZgwMYaeFTR4D6TfnvhgUnanXTzuz4WzmBQaG6MsMDEcueOpxnF5iwHZ+kfHaUypn
+                n5xefpbrzEYjZ3MGhiogDr/YYbxbjYHhrH6lYcY55zNgZzGcNWBgUL0Uctr3zLzTt08xMOScZmCYdnbl
+                qQMnpcxB8konQSQACVZG3A==''')
+        with self.subTest('navier-stokes-pressure'):
+            self.assertAlmostEqual64(args1['p'], '''
+                eNrbqjVZs1/ry/n48z1nSrW9L83RkTmneNZMO/TCOUNbMwDktQ3z''')
+
+    def test_compatible(self):
+        args0, args1 = main(nelems=3, etype='square', degree=2, reynolds=100, compatible=True, strongbc=False)
+        with self.subTest('stokes-velocity'):
+            self.assertAlmostEqual64(args0['u'], '''
+                eNpjYIAAvwvdBr9O2Zk90E8+rXQ6yxzGZ4CDTfr3z0H45hc2mjSagFgn9f1P15+G6Fc0PHSSgQEAx7kX
+                6A==''')
+        with self.subTest('stokes-pressure'):
+            self.assertAlmostEqual64(args0['p'], '''
+                eNoL1u+7NOfUR929ugvORxlU6W7V1TcUuyiif/PKCf1yUwDfRw2t''')
+        with self.subTest('navier-stokes-velocity'):
+            self.assertAlmostEqual64(args1['u'], '''
+                eNpjYICA1HNRRkGnZ5r26m86bX3awvyBftS5C6dOmDHAwWxDmbMzTUEsrfMrTA6YgFjKV53OOJ0FsR7o
+                F561OMnAAAC5tRfX''')
+        with self.subTest('navier-stokes-pressure'):
+            self.assertAlmostEqual64(args1['p'], '''
+                eNoz1VYx3HT6t16w/uKz73Uv6R7RNzx35swh7XdXrQ0TzADbMQ6l''')
+
+    def test_strong(self):
+        args0, args1 = main(nelems=3, etype='square', degree=2, reynolds=100, compatible=False, strongbc=True)
         with self.subTest('stokes-velocity'):
             self.assertAlmostEqual64(args0['u'], '''
-                eNpjYEAFEy++uHzs/EYjEFv73Hm9T2eVDGFyMWdfGJ/TVTGxMvQyqjxbAlYTZM7AsNWIxwhEG5hs0wTR
-                1Wd5DDTOHrx05OzmczC9Cud8riSccbrqYJR2dfMZ07M+hkznLpuongepB+Eyk/TzIHVbL4QbrLqw9Qyy
-                m+aee31awWALXEzP+NIZkB6Y/QAD2Dbr''')
+                eNpjYMAPDl2wNEg9p2D8QeOQcafBJ9OTJ6abB5lD5E6eVb348oyVkfSZf8YFZ6RNZp6+ZwiTO3KGgUEP
+                iNedYmDoPc3AsMGEgQGh77beyzPHzkqfYTpTcObp6Zmnlc7B5A5dOH4+9dyDMx807M50GvCdOnki8wRM
+                DhcAAEYiNtQ=''')
         with self.subTest('stokes-pressure'):
             self.assertAlmostEqual64(args0['p'], '''
-                eNoBIADf/1zNa81jzWHNs8w3zV3MUs2HzZrNzc26zanNd82qzgAAR/MTmQ==''')
+                eNoBHgDh/3fRlNYxy7PR0NVKz1ktVi1E1HowsdGJ07Qt/9PINA6QEUk=''')
         with self.subTest('navier-stokes-velocity'):
             self.assertAlmostEqual64(args1['u'], '''
-                eNpjYEAFvRcfXLa8eNsQxP5xLkBv17ktBjC5iHOfjbr1XxotNbAy0j4nbQQSCzJnYLA0lNIH0XLGi3VB
-                dNK5HoNjZ2frXj+77BxM775zfbq6Z+cYxBpd1vc543jursGSC0omz86D1IPwBJOzYDszL5oaLbgQcQbZ
-                TTZnec9svix4FsZXNbl9GqQHZj8AAcY1/g==''')
+                eNpjYMAPAs4H6M8zaDJOO91vKma628TihKx5kDlEbv2Z5fqFZ24aKZ/mNplmmGJy5eRbY5jcpdOGF+tP
+                /9BPPW1gXH6W2eSpkds5mBz72fvnUs/EnHt+etXpCWdZzgSd3W8Ck1M9L3zGGyi/4Pz80+ZnNp44c9L8
+                BEwOFwAA4RM3QA==''')
         with self.subTest('navier-stokes-pressure'):
             self.assertAlmostEqual64(args1['p'], '''
-                eNoBIADf/wXMDswMzAfMzMvry73LAcwIzCDM/ssJzDTM9MvRzAAAHqQRsw==''')
+                eNoBHgDh//ot489SzMEsntHezWTPAC+jL+XN/8wF02UxTc1JNhf0ELc=''')
diff --git a/nutils/export.py b/nutils/export.py
index c1494e958..d6d1b795e 100644
--- a/nutils/export.py
+++ b/nutils/export.py
@@ -47,34 +47,59 @@ def plotlines_(ax, xy, lines, **kwargs):
     return lc
 
 
-def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None):
-    if (tri is None) != (values is None):
+def _triplot_1d(ax, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None):
+    if plabel:
+        ax.set_xlabel(plabel)
+    if vlabel:
+        ax.set_ylabel(vlabel)
+    if hull is not None:
+        for x in points[hull[:, 0], 0]:
+            ax.axvline(x, color=linecolor, linewidth=linewidth)
+    if tri is not None:
+        plotlines_(ax, [points[:, 0], values], tri)
+    ax.autoscale(enable=True, axis='x', tight=True)
+    if clim is None:
+        ax.autoscale(enable=True, axis='y', tight=False)
+    else:
+        ax.set_ylim(clim)
+
+
+def _triplot_2d(ax, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None):
+    ax.set_aspect('equal')
+    if plabel:
+        ax.set_xlabel(plabel)
+        ax.set_ylabel(plabel)
+    if tri is not None:
+        im = ax.tripcolor(*points.T, tri, values, shading='gouraud', cmap=cmap, rasterized=True)
+        if clim is not None:
+            im.set_clim(clim)
+    else:
+        im = None
+    if hull is not None:
+        plotlines_(ax, points.T, hull, colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5)
+    ax.autoscale(enable=True, axis='both', tight=True)
+    return im
+
+
+def triplot(name, points, values=None, **kwargs):
+
+    if points.shape[1] == 1:
+        _triplot = _triplot_1d
+    elif points.shape[1] == 2:
+        _triplot = _triplot_2d
+    else:
+        raise Exception(f'invalid spatial dimension: {nd}')
+
+    if (kwargs.get('tri') is None) != (values is None):
         raise Exception('tri and values can only be specified jointly')
+
+    if not isinstance(name, str):
+        return _triplot(name, points, values, **kwargs)
+
     with mplfigure(name) as fig:
-        if points.shape[1] == 1:
-            ax = fig.add_subplot(111, xlabel=plabel, ylabel=vlabel)
-            if tri is not None:
-                plotlines_(ax, [points[:, 0], values], tri)
-            if hull is not None:
-                for x in points[hull[:, 0], 0]:
-                    ax.axvline(x, color=linecolor, linewidth=linewidth)
-            ax.autoscale(enable=True, axis='x', tight=True)
-            if clim is None:
-                ax.autoscale(enable=True, axis='y', tight=False)
-            else:
-                ax.set_ylim(clim)
-        elif points.shape[1] == 2:
-            ax = fig.add_subplot(111, xlabel=plabel, ylabel=plabel, aspect='equal')
-            if tri is not None:
-                im = ax.tripcolor(*points.T, tri, values, shading='gouraud', cmap=cmap, rasterized=True)
-                if clim is not None:
-                    im.set_clim(clim)
-                fig.colorbar(im, label=vlabel)
-            if hull is not None:
-                plotlines_(ax, points.T, hull, colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5)
-            ax.autoscale(enable=True, axis='both', tight=True)
-        else:
-            raise Exception('invalid spatial dimension: {}'.format(points.shape[1]))
+        im = _triplot(fig.add_subplot(111), points, values, **kwargs)
+        if im:
+            fig.colorbar(im, label=kwargs.get('vlabel'))
 
 
 @util.positional_only