diff --git a/examples/adaptivity.py b/examples/adaptivity.py
index 71d468e1c..b03d2caa5 100644
--- a/examples/adaptivity.py
+++ b/examples/adaptivity.py
@@ -53,18 +53,18 @@ def main(etype: str = 'square',
 
         if irefine:
             refdom = domain.refined
-            refbasis = refdom.basis(btype, degree=degree)
-            ns.add_field('vref', refbasis)
-            res = refdom.integral('∇_k(vref) ∇_k(u) dV' @ ns, degree=degree*2)
-            res -= refdom.boundary.integral('vref ∇_k(u) n_k dS' @ ns, degree=degree*2)
-            indicator = res.derivative('vref').eval(**args)
-            supp = refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
-            domain = domain.refined_by(refdom.transforms[supp])
+            ns.refbasis = refdom.basis(btype=btype, degree=degree)
+            res = refdom.integral('∇_k(refbasis_n) ∇_k(u) dV' @ ns, degree=degree*2)
+            res -= refdom.boundary.integral('refbasis_n ∇_k(u) n_k dS' @ ns, degree=degree*2)
+            indicator = numpy.square(res.eval(**args))
+            irefelems = ns.refbasis.get_support(indicator > indicator.mean())
+            domain = domain.refined_by(refdom.transforms[irefelems])
 
         ns = Namespace()
         ns.x = geom
         ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-        ns.add_field(('u', 'v'), domain.basis(btype, degree=degree))
+        ns.u = domain.field('u', btype=btype, degree=degree)
+        ns.v = domain.field('v', btype=btype, degree=degree)
         ns.uexact = exact
         ns.du = 'u - uexact'
 
diff --git a/examples/burgers.py b/examples/burgers.py
index 7e072efea..8f5839ff2 100644
--- a/examples/burgers.py
+++ b/examples/burgers.py
@@ -40,14 +40,16 @@ def main(nelems: int = 40,
     ns = Namespace()
     ns.x = geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-    ns.add_field(('u', 'u0', 'v'), domain.basis(btype, degree=degree))
-    ns.add_field(('t', 't0'))
-    ns.dudt = '(u - u0) / (t - t0)'
+    ns.u = domain.field('u', btype=btype, degree=degree)
+    ns.du = ns.u - function.replace_arguments(ns.u, 'u:u0')
+    ns.v = domain.field('v', btype=btype, degree=degree)
+    ns.t = function.field('t')
+    ns.dt = ns.t - function.field('t0')
     ns.f = '.5 u^2'
     ns.C = 1
     ns.uinit = 'exp(-25 x^2)'
 
-    res = domain.integral('(v dudt - ∇(v) f) dV' @ ns, degree=degree*2)
+    res = domain.integral('(v du / dt - ∇(v) f) dV' @ ns, degree=degree*2)
     res -= domain.interfaces.integral('[v] n ({f} - .5 C [u] n) dS' @ ns, degree=degree*2)
 
     sqr = domain.integral('(u - uinit)^2 dV' @ ns, degree=max(degree*2, 5))
diff --git a/examples/cahnhilliard.py b/examples/cahnhilliard.py
index 9485acb22..091f64c47 100644
--- a/examples/cahnhilliard.py
+++ b/examples/cahnhilliard.py
@@ -147,20 +147,19 @@ def main(size: Length = parse('10cm'),
         domain, geom = mesh.unitsquare(nelems, etype)
 
     ns = Namespace()
-    ns.x = size * geom
+    ns.x = geom * size
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-    ns.add_field(('φ', 'φ0'), domain.basis('std', degree=degree))
-    ns.add_field('η', domain.basis('std', degree=degree) * stens / epsilon) # basis scaling to give η the required unit
-    ns.add_field('dt')
-    ns.dt *= timestep
+    ns.φ = domain.field('φ', btype='std', degree=degree)
+    ns.dφ = ns.φ - function.replace_arguments(ns.φ, 'φ:φ0')
+    ns.η = domain.field('η', btype='std', degree=degree) * (stens / epsilon)
+    ns.dt = function.field('dt') * timestep
     ns.ε = epsilon
     ns.σ = stens
     ns.σmean = (wtensp + wtensn) / 2
     ns.σdiff = (wtensp - wtensn) / 2
     ns.σwall = 'σmean + φ σdiff'
-    ns.dφ = 'φ - φ0'
     ns.ψ = '.25 (φ^2 - 1)^2'
-    ns.δψ = '.25 dφ^2 (1 - φ0^2 / 6 - φ φ0 / 3 - φ^2 / 2)' if stable else '0'
+    ns.δψ = '.25 dφ^2 (1 - φ^2 + 2 φ dφ / 3 - dφ^2 / 6)' if stable else '0'
     ns.M = mobility
     ns.J_i = '-M ∇_i(η)'
     ns.v_i = 'φ J_i'
@@ -181,7 +180,7 @@ def main(size: Length = parse('10cm'),
     system = System(nrg / tol, trial='φ,η')
 
     numpy.random.seed(seed)
-    args = dict(φ=numpy.random.normal(0, .5, system.argshapes['φ'])) # initial condition
+    args = dict(φ=numpy.random.normal(0, .5, function.arguments_for(nrg)['φ'].shape)) # initial condition
 
     with log.iter.fraction('timestep', range(round(endtime / timestep))) as steps:
         for istep in steps:
diff --git a/examples/coil.py b/examples/coil.py
index 07ae8ec94..6425c8ff2 100644
--- a/examples/coil.py
+++ b/examples/coil.py
@@ -105,7 +105,8 @@ def main(nelems: int = 50,
     X = RZ * REV
     ns.x = ns.rz @ ns.rot
     ns.define_for('x', gradient='∇', jacobians=('dV', 'dS'), curl='curl')
-    ns.add_field(('A', 'Atest'), RZ.basis('spline', degree=degree, removedofs=[[0, -1], [-1]])[:,numpy.newaxis] * ns.eθ, dtype=complex)
+    ns.A = RZ.field('A', btype='spline', degree=degree, removedofs=[[0, -1], [-1]], dtype=complex) * ns.eθ
+    ns.Atest = function.replace_arguments(ns.A, 'A:Atest')
     ns.B_i = 'curl_ij(A_j)'
     ns.E_i = '-j ω A_i'
     ns.Jind_i = 'σ E_i'
@@ -123,7 +124,8 @@ def main(nelems: int = 50,
     # on a basis.
 
     ns.Borig = ns.B
-    ns.add_field(('B', 'Btest'), RZ.basis('spline', degree=degree), ns.rot, dtype=complex)
+    ns.B = function.field('B', RZ.basis('spline', degree=degree), ns.rot, dtype=complex)
+    ns.Btest = function.replace_arguments(ns.B, 'B:Btest')
     res = REV.integral(RZ.integral('Btest_i (B_i - Borig_i) dV' @ ns, degree=2*degree), degree=0)
     args = System(res, trial='B', test='Btest').solve(arguments=args)
 
diff --git a/examples/cylinderflow.py b/examples/cylinderflow.py
index 3a1f7ace1..f2fdd0bd7 100644
--- a/examples/cylinderflow.py
+++ b/examples/cylinderflow.py
@@ -115,11 +115,15 @@ def main(nelems: int = 99,
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
     J = ns.x.grad(geom)
     detJ = numpy.linalg.det(J)
-    ns.add_field(('u', 'u0', 'v'), function.vectorize([
+    ns.u = function.field('u', function.vectorize([
         domain.basis('spline', degree=(degree, degree-1), removedofs=((0,), None)),
         domain.basis('spline', degree=(degree-1, degree))]) @ J.T / detJ)
-    ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1) / detJ)
-    ns.add_field(('t', 'dt'))
+    ns.p = domain.field('p', btype='spline', degree=degree-1) / detJ
+    ns.v = function.replace_arguments(ns.u, 'u:v')
+    ns.q = function.replace_arguments(ns.p, 'p:q')
+    ns.t = function.field('t')
+    ns.du = ns.u - function.replace_arguments(ns.u, 'u:u0')
+    ns.dt = function.field('dt')
     ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
     ns.ω = 'ε_ij ∇_j(u_i)'
     ns.N = 10 * degree / elemangle  # Nitsche constant based on element size = elemangle/2
@@ -133,7 +137,7 @@ def main(nelems: int = 99,
     sqr = domain.integral('(.5 Σ_i (u_i - uinf_i)^2 - ∇_k(u_k) p) dV' @ ns, degree=degree*2)
     args = System(sqr, trial='u,p').solve(constrain=cons) # set initial condition to potential flow
 
-    res = domain.integral('v_i (u_i - u0_i) dV' @ ns, degree=degree*3)
+    res = domain.integral('v_i du_i dV' @ ns, degree=degree*3)
     res += domain.integral('(v_i ∇_j(u_i) u_j + ∇_j(v_i) σ_ij + q ∇_k(u_k)) dt dV' @ ns, degree=degree*3)
     res += domain.boundary['inner'].integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dt dS' @ ns, degree=degree*2)
     div = numpy.sqrt(abs(function.factor(domain.integral('∇_k(u_k)^2 dV' @ ns, degree=2)))) # L2 norm of velocity divergence
diff --git a/examples/drivencavity.py b/examples/drivencavity.py
index c110924af..cf70d42d2 100644
--- a/examples/drivencavity.py
+++ b/examples/drivencavity.py
@@ -100,15 +100,15 @@ def main(nelems: int = 32,
     ns.x = geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
     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:])
+        ns.u = domain.field('u', btype='std', degree=degree, shape=[2])
+        ns.p = domain.field('p', btype='std', degree=degree-1)
+        ns.ψ = domain.field('ψ', btype='std', degree=2)
     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.u = function.field('u', function.vectorize([domain.basis('spline', degree=p) for p in degree - 1 + numpy.eye(2, dtype=int)]))
+        ns.p = domain.field('p', btype='spline', degree=degree-1)
+        ns.ψ = domain.field('ψ', btype='spline', degree=degree)
+    ns.v = function.replace_arguments(ns.u, 'u:v')
+    ns.q = function.replace_arguments(ns.p, 'p:q')
     ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij'
 
     # weak formulation for Stokes flow, over-integrating for improved
@@ -119,6 +119,8 @@ def main(nelems: int = 32,
     # strong enforcement of non-penetrating boundary conditions
     sqr = domain.boundary.integral('(u_k n_k)^2 dS' @ ns, degree=degree*2)
     cons = System(sqr, trial='u').solve_constraints(droptol=1e-15)
+    cons['p'] = numpy.zeros(function.arguments_for(res)['p'].shape, dtype=bool)
+    cons['p'].flat[0] = True # point constraint
 
     if strongbc:
         # strong enforcement of tangential boundary conditions
@@ -156,7 +158,9 @@ def postprocess(domain, ns, **arguments):
 
     # reconstruct velocity streamlines
     sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(ψ))^2 dV' @ ns, degree=4)
-    arguments = System(sqr, trial='ψ').solve(arguments=arguments)
+    consψ = numpy.zeros(function.arguments_for(sqr)['ψ'].shape, dtype=bool)
+    consψ.flat[0] = True # point constraint
+    arguments = System(sqr, trial='ψ').solve(arguments=arguments, constrain={'ψ': consψ})
 
     bezier = domain.sample('bezier', 9)
     x, u, p, ψ = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'ψ'] @ ns, **arguments)
@@ -197,7 +201,7 @@ def test_baseline(self):
                 YcDGBdLMJR3Y/X+zdkhHHrEHM6lENt25+OU8OUi8PUn+klm87lacqiN4uQrZ4tUCLh3g4AFrV6Q/uctG
                 gQ==''')
         with self.subTest('stokes-pressure'):
-            self.assertAlmostEqual64(args0['p'], '''
+            self.assertAlmostEqual64(args0['p'][1:], '''
                 eNoBHgDh/+vVsNEXy6jTbiq1z7Av9C0mLJkw1NDTLEEtEC/xNAwED0s=''')
         with self.subTest('navier-stokes-velocity'):
             self.assertAlmostEqual64(args1['u'], '''
@@ -206,7 +210,7 @@ def test_baseline(self):
                 nM/5nf5xevqZxDOq5w4bCwLlOoD6XDV/n1t//s5ZvjPzTjmdDjx55+Slky/MGaDgHFB3vz4DgynQfS9O
                 A3WcBIkCAB7aSkk=''')
         with self.subTest('navier-stokes-pressure'):
-            self.assertAlmostEqual64(args1['p'], '''
+            self.assertAlmostEqual64(args1['p'][1:], '''
                 eNpz01W9oHVmuU7SJYtzgherdcr0n59dfiZT11yP97yCGQDN0Azu''')
 
     def test_mixed(self):
@@ -218,7 +222,7 @@ def test_mixed(self):
                 n+Y8k3RG+Mwio99nuoHqg4G48WzCmTignYUXDfXNzoedATuL4bMeA0Op9qczWqfXnTl2ioHhINAdHufv
                 ntx18qoZSH7FSRAJAB13Sc0=''')
         with self.subTest('stokes-pressure'):
-            self.assertAlmostEqual64(args0['p'], '''
+            self.assertAlmostEqual64(args0['p'][1:], '''
                 eNp7pKl+nf1KznmxS62ns/W+az/TNTL4ondU1/46t6GKKQDiJg1H''')
         with self.subTest('navier-stokes-velocity'):
             self.assertAlmostEqual64(args1['u'], '''
@@ -227,7 +231,7 @@ def test_mixed(self):
                 n5xefpbrzEYjZ3MGhiogDr/YYbxbjYHhrH6lYcY55zNgZzGcNWBgUL0Uctr3zLzTt08xMOScZmCYdnbl
                 qQMnpcxB8konQSQACVZG3A==''')
         with self.subTest('navier-stokes-pressure'):
-            self.assertAlmostEqual64(args1['p'], '''
+            self.assertAlmostEqual64(args1['p'][1:], '''
                 eNrbqjVZs1/ry/n48z1nSrW9L83RkTmneNZMO/TCOUNbMwDktQ3z''')
 
     def test_compatible(self):
@@ -237,14 +241,14 @@ def test_compatible(self):
                 eNpjYIAAvwvdBr9O2Zk90E8+rXQ6yxzGZ4CDTfr3z0H45hc2mjSagFgn9f1P15+G6Fc0PHSSgQEAx7kX
                 6A==''')
         with self.subTest('stokes-pressure'):
-            self.assertAlmostEqual64(args0['p'], '''
+            self.assertAlmostEqual64(args0['p'][1:], '''
                 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'], '''
+            self.assertAlmostEqual64(args1['p'][1:], '''
                 eNoz1VYx3HT6t16w/uKz73Uv6R7RNzx35swh7XdXrQ0TzADbMQ6l''')
 
     def test_strong(self):
@@ -255,7 +259,7 @@ def test_strong(self):
                 iNedYmDoPc3AsMGEgQGh77beyzPHzkqfYTpTcObp6Zmnlc7B5A5dOH4+9dyDMx807M50GvCdOnki8wRM
                 DhcAAEYiNtQ=''')
         with self.subTest('stokes-pressure'):
-            self.assertAlmostEqual64(args0['p'], '''
+            self.assertAlmostEqual64(args0['p'][1:], '''
                 eNoBHgDh/3fRlNYxy7PR0NVKz1ktVi1E1HowsdGJ07Qt/9PINA6QEUk=''')
         with self.subTest('navier-stokes-velocity'):
             self.assertAlmostEqual64(args1['u'], '''
@@ -263,7 +267,7 @@ def test_strong(self):
                 /9BPPW1gXH6W2eSpkds5mBz72fvnUs/EnHt+etXpCWdZzgSd3W8Ck1M9L3zGGyi/4Pz80+ZnNp44c9L8
                 BEwOFwAA4RM3QA==''')
         with self.subTest('navier-stokes-pressure'):
-            self.assertAlmostEqual64(args1['p'], '''
+            self.assertAlmostEqual64(args1['p'][1:], '''
                 eNoBHgDh//ot489SzMEsntHezWTPAC+jL+XN/8wF02UxTc1JNhf0ELc=''')
 
 
diff --git a/examples/elasticity.py b/examples/elasticity.py
index 7cc82879a..f08d9dfa2 100644
--- a/examples/elasticity.py
+++ b/examples/elasticity.py
@@ -41,7 +41,7 @@ def main(nelems: int = 24,
     ns.δ = function.eye(domain.ndims)
     ns.x = geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-    ns.add_field(('u', 't'), domain.basis(btype, degree=degree), shape=(2,))
+    ns.u = domain.field('u', btype=btype, degree=degree, shape=[2])
     ns.X_i = 'x_i + u_i'
     ns.λ = 1
     ns.μ = .5/poisson - 1
@@ -61,8 +61,9 @@ def main(nelems: int = 24,
     if direct:
         ns.t_i = 'σ_ij n_j' # <-- this is an inadmissible boundary term
     else:
-        energy -= domain.boundary['top'].integral('u_i t_i dS' @ ns, degree=degree*2)
-        args = System(energy, trial='t', test='u').solve(constrain={'t': numpy.isnan(cons['u'])}, arguments=args)
+        ns.t = domain.field('t', btype=btype, degree=degree, shape=[2])
+        system = System(energy - domain.boundary['top'].integral('u_i t_i dS' @ ns, degree=degree*2), trial='t', test='u')
+        args = system.solve(constrain={'t': numpy.isnan(cons['u'])}, arguments=args)
     F = domain.boundary['top'].integrate('t_i dS' @ ns, degree=degree*2, arguments=args)
     log.user('total clamping force:', F)
 
diff --git a/examples/finitestrain.py b/examples/finitestrain.py
index 058d5bd6b..5de553cab 100644
--- a/examples/finitestrain.py
+++ b/examples/finitestrain.py
@@ -51,7 +51,7 @@ def main(nelems: int = 20,
     ns.angle = angle * numpy.pi / 180
     ns.λ = 2 * poisson
     ns.μ = 1 - 2 * poisson
-    ns.add_field('u', domain.basis(btype, degree=degree), shape=(domain.ndims,))
+    ns.u = domain.field('u', btype=btype, degree=degree, shape=[domain.ndims])
     ns.x_i = 'X_i + u_i'
     ns.ε_ij = '.5 (∇_j(u_i) + ∇_i(u_j))'
     ns.energy = 'λ ε_ii ε_jj + 2 μ ε_ij ε_ij'
diff --git a/examples/laplace.py b/examples/laplace.py
index 0c52ddfbc..783a2cd5e 100644
--- a/examples/laplace.py
+++ b/examples/laplace.py
@@ -49,13 +49,13 @@ def main(nelems: int = 10,
 
     # To be able to write index based tensor contractions, we need to bundle
     # all relevant functions together in a namespace. Here we add the geometry
-    # `x`, a test function `v`, and the solution `u`. The latter two are formed
-    # by contracting a basis with function arguments of the same name.
+    # `x`, a test function `v`, and the solution `u`.
 
     ns = Namespace()
     ns.x = geom
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
-    ns.add_field(('u', 'v'), domain.basis(btype, degree=degree))
+    ns.u = domain.field('u', btype=btype, degree=degree)
+    ns.v = domain.field('v', btype=btype, degree=degree)
 
     # We are now ready to implement the Laplace equation. In weak form, the
     # solution is a scalar field `u` for which ∫_Ω ∇v·∇u - ∫_Γn v f = 0 ∀ v.
diff --git a/examples/platewithhole.py b/examples/platewithhole.py
index 5f8845858..ed86b0b32 100644
--- a/examples/platewithhole.py
+++ b/examples/platewithhole.py
@@ -75,7 +75,7 @@ def generate(self, radius):
         if self.nrefine:
             topo = topo.refine(self.nrefine)
             bsplinebasis = topo.basis('spline', degree=2)
-            sqr = topo.integral((function.dotarg('w', bsplinebasis) - weightfunc)**2, degree=9)
+            sqr = topo.integral((function.field('w', bsplinebasis) - weightfunc)**2, degree=9)
             controlweights = System(sqr, trial='w').solve()['w']
             nurbsbasis = bsplinebasis * controlweights / weightfunc
         return topo.withboundary(hole='left', sym='top,bottom', far='right'), geom, nurbsbasis, 5
@@ -121,7 +121,8 @@ def main(mode: Union[FCM, NURBS] = NURBS(),
     ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
     ns.λ = 2 * poisson
     ns.μ = 1 - poisson
-    ns.add_field(('u', 'v'), basis, shape=(2,))
+    ns.u = function.field('u', basis, shape=[2])
+    ns.v = function.field('v', basis, shape=[2])
     ns.X_i = 'x_i + u_i'
     ns.ε_ij = '(∇_j(u_i) + ∇_i(u_j)) / 2'
     ns.σ_ij = 'λ ε_kk δ_ij + 2 μ ε_ij'
diff --git a/examples/poisson.py b/examples/poisson.py
index 10081a1d0..8c5dc4e8e 100644
--- a/examples/poisson.py
+++ b/examples/poisson.py
@@ -20,7 +20,7 @@ def main(nelems: int = 32):
     '''
 
     topo, x = mesh.unitsquare(nelems, etype='square')
-    u = function.dotarg('u', topo.basis('std', degree=1))
+    u = topo.field('u', btype='std', degree=1)
     g = u.grad(x)
     J = function.J(x)
 
diff --git a/examples/torsion.py b/examples/torsion.py
index 9d0a84d9a..4826cc9f5 100644
--- a/examples/torsion.py
+++ b/examples/torsion.py
@@ -66,7 +66,7 @@ def main(length: float = 2*np.pi,
     zgrid = length * np.linspace(-.5, .5, round(length / elemsize)+1)
     θgrid = np.linspace(-np.pi, np.pi, round(2 * np.pi / elemsize)+1)
     cylinder, (z, θ) = mesh.rectilinear([zgrid, θgrid], periodic=(1,))
-    φ = θ - (z / length * np.pi / 180) * function.Argument('φ', shape=())
+    φ = θ - (z / length * np.pi / 180) * function.field('φ')
     if trim:
         cylinder = cylinder.trim(θ**2 + z**2 - trim**2, maxrefine=2)
     extrusion, r = mesh.line([1 - thickness/2, 1 + thickness/2], space='T')
@@ -77,8 +77,7 @@ def main(length: float = 2*np.pi,
     ns.X = np.stack([z, r * np.sin(θ), r * np.cos(θ)]) # reference geometry
     ns.Xφ = np.stack([z * stretch, r * np.sin(φ), r * np.cos(φ)])
     ns.define_for('X', gradient='∇', jacobians=('dV',))
-    ns.add_field('u', topo.basis('spline', degree=degree,
-        removedofs=((0,-1),None,None)), shape=(3,)) # deformation, clamped
+    ns.u = topo.field('u', btype='spline', degree=degree, removedofs=((0,-1),None,None), shape=[3]) # deformation, clamped
     ns.x_i = 'Xφ_i + u_i' # deformed geometry
     ns.F_ij = '∇_j(x_i)'
     ns.J = np.linalg.det(ns.F)
diff --git a/nutils/SI.py b/nutils/SI.py
index 284f7bdf7..43ec2cc0c 100644
--- a/nutils/SI.py
+++ b/nutils/SI.py
@@ -406,11 +406,16 @@ def __evaluate(op, *args, **kwargs):
         dims, args = zip(*Quantity.__unpack(*args))
         return tuple(dim.wrap(ret) for (dim, ret) in zip(dims, op(*args, **kwargs)))
 
-    @register('dotarg')
-    def __dotarg(op, *args, **kwargs):
+    @register('field')
+    def __field(op, *args, **kwargs):
         dims, args = zip(*Quantity.__unpack(*args)) # we abuse the fact that unpack str returns dimensionless
         return functools.reduce(operator.mul, dims).wrap(op(*args, **kwargs))
 
+    @register('arguments_for')
+    def __attribute(op, *args, **kwargs):
+        __dims, args = zip(*Quantity.__unpack(*args))
+        return op(*args, **kwargs)
+
     del register
 
     ## DEFINE OPERATORS
diff --git a/nutils/__init__.py b/nutils/__init__.py
index b2e00bb05..3ee2f104a 100644
--- a/nutils/__init__.py
+++ b/nutils/__init__.py
@@ -1,4 +1,4 @@
 'Numerical Utilities for Finite Element Analysis'
 
-__version__ = version = '9a46'
+__version__ = version = '9a47'
 version_name = 'jook-sing'
diff --git a/nutils/evaluable.py b/nutils/evaluable.py
index b33de5f3c..3ec07028c 100644
--- a/nutils/evaluable.py
+++ b/nutils/evaluable.py
@@ -3394,27 +3394,40 @@ def evalf(func, *args):
         return numeric.accumulate(func, reshaped_indices, shape)
 
     def _compile_with_out(self, builder, out, out_block_id, mode):
+        # Compiles to an assignment (or in place addition) of the form:
+        #
+        #     out[index1, index2, ...] = func.transpose(trans)
+        #
+        # The order of the indices is unchanged, but ranges are converted to
+        # slices and any other indices (dubbed 'advanced') reshaped in similar
+        # fashion to numpy._ix to index the cross product. The right hand side
+        # is transposed if a slice operation separates two advanced indices to
+        # match Numpy's advanced indexing rules.
         assert mode in ('iadd', 'assign')
         if mode == 'assign':
             builder.get_block_for_evaluable(self, block_id=out_block_id, comment='zero').array_fill_zeros(out)
-        compiled_indices = []
         advanced_ndim = builtins.sum(index.ndim for index in self.indices if not isinstance(index, Range))
+        compiled_indices = []
+        trans = [] # axes of func corresponding to advanced indices
         i = 0
         for index in self.indices:
+            j = i + index.ndim
             if isinstance(index, Range):
                 n = builder.compile(index.shape[0])
-                compiled_indices.append(_pyast.Variable('slice').call(n))
+                compiled_index = _pyast.Variable('slice').call(n)
             else:
-                j = i + index.ndim
+                prefix = len(trans)
+                trans.extend(range(i, j))
+                suffix = advanced_ndim - len(trans)
                 compiled_index = builder.compile(index)
                 if index.ndim < advanced_ndim:
-                    compiled_index = compiled_index.get_item(_pyast.Raw(','.join(['None'] * i + [':'] * index.ndim + ['None'] * (advanced_ndim-j))))
-                compiled_indices.append(compiled_index)
-                i = j
-        assert i == advanced_ndim
+                    compiled_index = compiled_index.get_item(_pyast.Raw(','.join(['None'] * prefix + [':'] * index.ndim + ['None'] * suffix)))
+            compiled_indices.append(compiled_index)
+            i = j
+        assert i == self.func.ndim
+        assert len(trans) == advanced_ndim
         compiled_func = builder.compile(self.func)
-        trans = [i for i, index in enumerate(self.indices) if not isinstance(index, Range)]
-        if advanced_ndim and any(numpy.diff(trans) > 1):
+        if advanced_ndim > 1 and trans[-1] - trans[0] != advanced_ndim - 1: # trans is noncontiguous
             # see https://numpy.org/doc/stable/user/basics.indexing.html#combining-advanced-and-basic-indexing
             trans.extend(i for i, index in enumerate(self.indices) if isinstance(index, Range))
             compiled_func = compiled_func.get_attr('transpose').call(*[_pyast.LiteralInt(i) for i in trans])
@@ -3422,22 +3435,27 @@ def _compile_with_out(self, builder, out, out_block_id, mode):
 
     def _optimized_for_numpy(self):
         if isinstance(self.func, Assemble):
-            nontrivial_indices = []
-            for f in self.func, self: # order is important to maintain insertion order of ndim=0 indices
-                offset = 0
-                for index in f.indices:
-                    if index != Range(f.shape[offset]):
-                        nontrivial_indices.append((offset, index))
-                    offset += index.ndim
-                assert offset == f.func.ndim
-            nontrivial_indices.sort(key=lambda item: item[0]) # order of equal offsets is maintained
-            indices = []
-            for i, index in nontrivial_indices:
-                if i < len(indices):
-                    return # inflations overlap
-                while i > len(indices):
-                    indices.append(Range(self.shape[len(indices)]))
-                indices.append(index)
+            indices = list(self.indices)
+            # We aim to merge the indices from the nested Assemble operations
+            # if they are separable, i.e. preceded by or following on full
+            # slices, by replacing Range instances in indices by the
+            # corresponding index from self.func.
+            for i, index in enumerate(self.func.indices):
+                if index != Range(self.func.shape[i]): # non-trivial index of self.func
+                    # we need to account for the different axis numberings
+                    # between self and self.func to find the right insertion
+                    # point.
+                    ax1 = 0 # axis of self.func
+                    ax2 = 0 # axis of self
+                    while ax1 < i: # find ax1, ax2 corresponding to i
+                        ax1 += indices[ax2].ndim
+                        ax2 += 1
+                    if ax1 != i or ax2 >= self.ndim or indices[ax2] != Range(self.shape[ax2]):
+                        # Any nontrivial nesting scenario would have been
+                        # handled by Inflate if possible, so we simply bail out
+                        # at the first sign of difficulty.
+                        return
+                    indices[ax2] = index # merge!
             return Assemble(self.func.func, tuple(indices), self.shape)
 
 
diff --git a/nutils/expression_v2.py b/nutils/expression_v2.py
index ade879738..be90be1d8 100644
--- a/nutils/expression_v2.py
+++ b/nutils/expression_v2.py
@@ -695,8 +695,8 @@ def define_for(self, __name: str, *, gradient: Optional[str] = None, curl: Optio
         >>> topo, ns.x = mesh.rectilinear([2, 2])
         >>> ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
         >>> ns.basis = topo.basis('spline', degree=1)
-        >>> ns.u = function.dotarg('u', ns.basis)
-        >>> ns.v = function.dotarg('v', ns.basis)
+        >>> ns.u = function.field('u', ns.basis)
+        >>> ns.v = function.field('v', ns.basis)
         >>> res = topo.integral('-∇_i(v) ∇_i(u) dV' @ ns, degree=2)
         >>> res += topo.boundary.integral('∇_i(v) u n_i dS' @ ns, degree=2)
         '''
@@ -718,7 +718,7 @@ def define_for(self, __name: str, *, gradient: Optional[str] = None, curl: Optio
             setattr(self, jacobian, function.jacobian(geom, numpy.size(geom) - i))
 
     def add_field(self, __names: Union[str, Sequence[str]], *__bases, shape: Tuple[int, ...] = (), dtype: function.DType = float):
-        '''Add field(s) of the form ns.u = function.dotarg('u', ...)
+        '''Add field(s) of the form ns.u = function.field('u', ...)
 
         Parameters
         ----------
@@ -733,7 +733,7 @@ def add_field(self, __names: Union[str, Sequence[str]], *__bases, shape: Tuple[i
         '''
 
         for name in (__names,) if isinstance(__names, str) else __names:
-            setattr(self, name, function.dotarg(name, *__bases, shape=shape, dtype=dtype))
+            setattr(self, name, function.field(name, *__bases, shape=shape, dtype=dtype))
 
     def copy_(self, **replacements: Mapping[str, function.Array]) -> 'Namespace':
         '''Return a copy of this namespace.
diff --git a/nutils/function.py b/nutils/function.py
index 3f74942a0..e4b8b9a0a 100644
--- a/nutils/function.py
+++ b/nutils/function.py
@@ -224,7 +224,7 @@ def __init__(self, shape: Shape, dtype: DType, spaces: FrozenSet[str], arguments
         self.shape = tuple(sh.__index__() for sh in shape)
         self.dtype = dtype
         self.spaces = frozenset(spaces)
-        self.arguments = dict(arguments)
+        self._arguments = dict(arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         raise NotImplementedError
@@ -236,8 +236,8 @@ def as_evaluable_array(self) -> evaluable.Array:
         return self.lower(LowerArgs((), {}, {}))
 
     def __index__(self):
-        if self.arguments or self.spaces:
-            raise ValueError('cannot convert non-constant array to index: arguments={}'.format(','.join(self.arguments)))
+        if self._arguments or self.spaces:
+            raise ValueError('cannot convert non-constant array to index: arguments={}'.format(','.join(self._arguments)))
         elif self.ndim:
             raise ValueError('cannot convert non-scalar array to index: shape={}'.format(self.shape))
         elif self.dtype != int:
@@ -483,11 +483,17 @@ def replace(self, __arguments: Mapping[str, IntoArray]) -> 'Array':
 
     def contains(self, __name: str) -> bool:
         'Test if target occurs in this function.'
-        return __name in self.arguments
+        return __name in self._arguments
+
+    @property
+    def arguments(self) -> Mapping[str, Tuple[Shape, DType]]:
+        warnings.deprecation('array.arguments is deprecated and will be removed in Nutils 10, please use function.arguments_for(array) instead')
+        return self._arguments.copy()
 
     @property
     def argshapes(self) -> Mapping[str, Tuple[int, ...]]:
-        return {name: shape for name, (shape, dtype) in self.arguments.items()}
+        warnings.deprecation("array.argshapes[...] is deprecated and will be removed in Nutils 10, please use function.arguments_for(array)[...].shape instead")
+        return {name: shape for name, (shape, dtype) in self._arguments.items()}
 
     def conjugate(self):
         '''Return the complex conjugate, elementwise.
@@ -684,7 +690,7 @@ def __init__(self, args: Iterable[Any], shape: Tuple[int], dtype: DType, npointw
         self._args = args
         self._npointwise = npointwise
         spaces = frozenset(space for arg in args if isinstance(arg, Array) for space in arg.spaces)
-        arguments = _join_arguments(arg.arguments for arg in args if isinstance(arg, Array))
+        arguments = _join_arguments(arg._arguments for arg in args if isinstance(arg, Array))
         super().__init__(shape=(*points_shape, *shape), dtype=dtype, spaces=spaces, arguments=arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
@@ -692,7 +698,7 @@ def lower(self, args: LowerArgs) -> evaluable.Array:
         add_points_shape = tuple(map(evaluable.asarray, self.shape[:self._npointwise]))
         points_shape = args.points_shape + add_points_shape
         coordinates = {space: evaluable.Transpose.to_end(evaluable.appendaxes(coords, add_points_shape), coords.ndim-1) for space, coords in args.coordinates.items()}
-        return _CustomEvaluable(type(self).__name__, self.evalf, self.partial_derivative, evalargs, self.shape[self._npointwise:], self.dtype, self.spaces, types.frozendict(self.arguments), LowerArgs(points_shape, types.frozendict(args.transform_chains), types.frozendict(coordinates)))
+        return _CustomEvaluable(type(self).__name__, self.evalf, self.partial_derivative, evalargs, self.shape[self._npointwise:], self.dtype, self.spaces, types.frozendict(self._arguments), LowerArgs(points_shape, types.frozendict(args.transform_chains), types.frozendict(coordinates)))
 
     @classmethod
     def evalf(cls, *args: Any) -> numpy.ndarray:
@@ -833,7 +839,7 @@ class _WithoutPoints:
     def __init__(self, __arg: Array) -> None:
         self._arg = __arg
         self.spaces = __arg.spaces
-        self.arguments = __arg.arguments
+        self._arguments = __arg._arguments
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         return self._arg.lower(LowerArgs((), args.transform_chains, {}))
@@ -851,7 +857,7 @@ def __init__(self, lower: Callable[..., evaluable.Array], *args: Lowerable, shap
         self._args = args
         assert all(hasattr(arg, 'lower') for arg in self._args)
         spaces = frozenset(space for arg in args for space in arg.spaces)
-        arguments = _join_arguments(arg.arguments for arg in self._args)
+        arguments = _join_arguments(arg._arguments for arg in self._args)
         super().__init__(shape, dtype, spaces, arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
@@ -911,25 +917,13 @@ class _Replace(Array):
     def __init__(self, arg: Array, replacements: Dict[str, Array]) -> None:
         self._arg = arg
         self._replacements = {}
+        for old, new in _argument_to_array(replacements, arg):
+            if new.spaces:
+                raise ValueError(f'replacement functions cannot be bound to a space, but replacement for Argument {old.name!r} is bound to {", ".join(new.spaces)}.')
+            self._replacements[old.name] = new
         # Build arguments map with replacements.
-        unreplaced = {}
-        for name, (shape, dtype) in arg.arguments.items():
-            if name in replacements:
-                replacement = replacements[name]
-                if isinstance(replacement, str):
-                    replacement = Argument(replacement, shape, dtype)
-                else:
-                    replacement = Array.cast(replacement)
-                    if replacement.shape != shape:
-                        raise ValueError(f'Argument {name!r} has shape {shape} but the replacement has shape {replacement.shape}.')
-                    if replacement.dtype != dtype:
-                        raise ValueError(f'Argument {name!r} has dtype {dtype.__name__} but the replacement has dtype {replacement.dtype.__name__}.')
-                    if replacement.spaces:
-                        raise ValueError(f'replacement functions cannot contain spaces, but replacement for Argument {name!r} contains space {", ".join(replacement.spaces)}.')
-                self._replacements[name] = replacement
-            else:
-                unreplaced[name] = shape, dtype
-        arguments = _join_arguments([unreplaced] + [replacement.arguments for replacement in self._replacements.values()])
+        unreplaced = {name: shape_dtype for name, shape_dtype in arg._arguments.items() if name not in replacements}
+        arguments = _join_arguments([unreplaced] + [replacement._arguments for replacement in self._replacements.values()])
         super().__init__(arg.shape, arg.dtype, arg.spaces, arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
@@ -962,7 +956,7 @@ def to_end(cls, array: Array, *axes: int) -> Array:
     def __init__(self, arg: Array, axes: Tuple[int, ...]) -> None:
         self._arg = arg
         self._axes = tuple(n.__index__() for n in axes)
-        super().__init__(tuple(arg.shape[axis] for axis in axes), arg.dtype, arg.spaces, arg.arguments)
+        super().__init__(tuple(arg.shape[axis] for axis in axes), arg.dtype, arg.spaces, arg._arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         arg = self._arg.lower(args)
@@ -976,7 +970,7 @@ class _Opposite(Array):
     def __init__(self, arg: Array, space: str) -> None:
         self._arg = arg
         self._space = space
-        super().__init__(arg.shape, arg.dtype, arg.spaces, arg.arguments)
+        super().__init__(arg.shape, arg.dtype, arg.spaces, arg._arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         oppargs = LowerArgs(args.points_shape, dict(args.transform_chains), args.coordinates)
@@ -1043,7 +1037,7 @@ def __init__(self, arg: Array, var: Argument) -> None:
         self._arg = arg
         self._var = var
         self._eval_var = evaluable.Argument(var.name, tuple(evaluable.constant(n) for n in var.shape), var.dtype)
-        arguments = _join_arguments((arg.arguments, var.arguments))
+        arguments = _join_arguments((arg._arguments, var._arguments))
         super().__init__(arg.shape+var.shape, complex if var.dtype == complex else arg.dtype, arg.spaces | var.spaces, arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
@@ -1068,7 +1062,7 @@ def __init__(self, func: Array, geom: Array) -> None:
         common_shape = broadcast_shapes(func.shape, geom.shape[:-1])
         self._func = numpy.broadcast_to(func, common_shape)
         self._geom = numpy.broadcast_to(geom, (*common_shape, geom.shape[-1]))
-        arguments = _join_arguments((func.arguments, geom.arguments))
+        arguments = _join_arguments((func._arguments, geom._arguments))
         super().__init__(self._geom.shape, complex if func.dtype == complex else float, func.spaces | geom.spaces, arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
@@ -1094,7 +1088,7 @@ def __init__(self, func: Array, geom: Array) -> None:
         common_shape = broadcast_shapes(func.shape, geom.shape[:-1])
         self._func = numpy.broadcast_to(func, common_shape)
         self._geom = numpy.broadcast_to(geom, (*common_shape, geom.shape[-1]))
-        arguments = _join_arguments((func.arguments, geom.arguments))
+        arguments = _join_arguments((func._arguments, geom._arguments))
         super().__init__(self._geom.shape, complex if func.dtype == complex else float, func.spaces | geom.spaces, arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
@@ -1124,7 +1118,7 @@ def __init__(self, geom: Array, tip_dim: Optional[int] = None) -> None:
                              'not greater than the dimension of the geometry.')
         self._tip_dim = tip_dim
         self._geom = geom
-        super().__init__((), float, geom.spaces, geom.arguments)
+        super().__init__((), float, geom.spaces, geom._arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         geom = self._geom.lower(args)
@@ -1145,7 +1139,7 @@ class _Normal(Array):
     def __init__(self, geom: Array) -> None:
         self._geom = geom
         assert geom.dtype == float
-        super().__init__(geom.shape, float, geom.spaces, geom.arguments)
+        super().__init__(geom.shape, float, geom.spaces, geom._arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         geom = self._geom.lower(args)
@@ -1180,7 +1174,7 @@ class _ExteriorNormal(Array):
     def __init__(self, rgrad: Array) -> None:
         assert rgrad.dtype == float and rgrad.shape[-2] == rgrad.shape[-1] + 1
         self._rgrad = rgrad
-        super().__init__(rgrad.shape[:-1], float, rgrad.spaces, rgrad.arguments)
+        super().__init__(rgrad.shape[:-1], float, rgrad.spaces, rgrad._arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         rgrad = self._rgrad.lower(args)
@@ -1207,7 +1201,7 @@ def __init__(self, __arrays: Sequence[IntoArray], axis: int) -> None:
             shape=(*shape0[:self.axis], builtins.sum(array.shape[self.axis] for array in self.arrays), *shape0[self.axis+1:]),
             dtype=self.arrays[0].dtype,
             spaces=functools.reduce(operator.or_, (array.spaces for array in self.arrays)),
-            arguments=_join_arguments(array.arguments for array in self.arrays))
+            arguments=_join_arguments(array._arguments for array in self.arrays))
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         return util.sum(evaluable._inflate(array.lower(args), evaluable.Range(evaluable.constant(array.shape[self.axis])) + offset, evaluable.constant(self.shape[self.axis]), self.axis-self.ndim)
@@ -1687,6 +1681,53 @@ def kronecker(__array: IntoArray, axis: int, length: int, pos: IntoArray) -> Arr
     return _Transpose.from_end(scatter(__array, length, pos), axis)
 
 
+def _argument_to_array(d: Any, array: Array) -> Iterable[Tuple[Argument, Array]]:
+    '''Helper function for argument replacement.
+
+    Given dictionary-like input, along with an array to look up argument names,
+    yield ``(arg, new)`` tuples where ``arg`` is an ``Argument`` and ``new`` an
+    ``Array`` of equal shape and dtype. A ``ValueError`` is raised if any key
+    cannot be made a valid argument for the provided array or any value cannot
+    be made an array of the same shape and dtype.
+
+    Dictionary-like input is either an actual dictionary, a sequence of (key,
+    value) tuples or strings, or a string, where any string is replaced by the
+    relevant Argument object. The following inputs are all equivalent:
+
+        d = {'u': 'v', 'p': 'q'}
+        d = 'u:v,p:q'
+        d = ('u:v', 'p:q')
+        d = [('u', 'v'), ('p', 'q')]
+        d = {'u': Argument('v', ...), 'p': Argument('q', ...)}
+        d = [(Argument('u', ...), Argument('v', ...)), 'p:q']
+    '''
+
+    for item in d.split(',') if isinstance(d, str) else d.items() if isinstance(d, dict) else d:
+        arg, new = item.split(':', 1) if isinstance(item, str) else item
+
+        if isinstance(arg, str):
+            if arg not in array._arguments:
+                continue
+            arg = Argument(arg, *array._arguments[arg])
+        elif not isinstance(arg, Argument):
+            raise ValueError('Key must be string or argument')
+        elif arg.name not in arguments:
+            continue
+        elif array._arguments[arg.name] != (arg.shape, arg.dtype):
+            raise ValueError(f'Argument {arg.name!r} has wrong shape or dtype')
+
+        if isinstance(new, str):
+            new = Argument(new, arg.shape, arg.dtype)
+        else:
+            new = Array.cast(new)
+            if new.shape != arg.shape:
+                raise ValueError(f'Argument {arg.name!r} has shape {arg.shape} but the replacement has shape {new.shape}.')
+            elif new.dtype != arg.dtype:
+                raise ValueError(f'Argument {arg.name!r} has dtype {arg.dtype.__name__} but the replacement has dtype {new.dtype.__name__}.')
+
+        yield arg, new
+
+
 @nutils_dispatch
 def replace_arguments(__array: IntoArray, __arguments: Mapping[str, Union[IntoArray, str]]) -> Array:
     '''Replace arguments with :class:`Array` objects.
@@ -1735,15 +1776,8 @@ def linearize(__array: IntoArray, __arguments: Union[str, Dict[str, str], Iterab
     '''
 
     array = Array.cast(__array)
-    args = __arguments.split(',') if isinstance(__arguments, str) \
-      else __arguments.items() if isinstance(__arguments, dict) \
-      else __arguments
-    parts = []
-    for kv in args:
-        k, v = kv.split(':', 1) if isinstance(kv, str) else kv
-        f = derivative(array, k)
-        parts.append(numpy.sum(f * Argument(v, *array.arguments[k]), tuple(range(array.ndim, f.ndim))))
-    return util.sum(parts)
+    return util.sum(numpy.sum(derivative(array, arg) * lin, array.ndim + numpy.arange(arg.ndim))
+        for arg, lin in _argument_to_array(__arguments, array))
 
 
 def broadcast_arrays(*arrays: IntoArray) -> Tuple[Array, ...]:
@@ -1824,14 +1858,14 @@ def derivative(__arg: IntoArray, __var: Union[str, 'Argument']) -> Array:
 
     arg = Array.cast(__arg)
     if isinstance(__var, str):
-        if __var not in arg.arguments:
+        if __var not in arg._arguments:
             raise ValueError('no such argument: {}'.format(__var))
-        shape, dtype = arg.arguments[__var]
+        shape, dtype = arg._arguments[__var]
         __var = Argument(__var, shape, dtype=dtype)
     elif not isinstance(__var, Argument):
         raise ValueError('Expected an instance of `Argument` as second argument of `derivative` but got a `{}.{}`.'.format(type(__var).__module__, type(__var).__qualname__))
-    if __var.name in arg.arguments:
-        shape, dtype = arg.arguments[__var.name]
+    if __var.name in arg._arguments:
+        shape, dtype = arg._arguments[__var.name]
         if __var.shape != shape:
             raise ValueError('Argument {!r} has shape {} in the function, but the derivative to {!r} with shape {} was requested.'.format(__var.name, shape, __var.name, __var.shape))
         if __var.dtype != dtype:
@@ -2318,6 +2352,7 @@ def vectorize(args: Sequence[IntoArray]) -> Array:
     :class:`Array`
     '''
 
+    args = tuple(args)
     return numpy.concatenate([kronecker(arg, axis=-1, length=len(args), pos=iarg) for iarg, arg in enumerate(args)])
 
 
@@ -2339,8 +2374,14 @@ def rotmat(__arg: IntoArray) -> Array:
     return Array.cast(numpy.stack([trignormal(__arg), trigtangent(__arg)], 0))
 
 
+def dotarg(*args, **kwargs):
+    '''Alias for :func:`field`.'''
+
+    return field(*args, **kwargs)
+
+
 @nutils_dispatch
-def dotarg(__argname: str, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtype: DType = float) -> Array:
+def field(name: str, /, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtype: DType = float) -> Array:
     '''Return the inner product of the first axes of the given arrays with an argument with the given name.
 
     An argument with shape ``(arrays[0].shape[0], ..., arrays[-1].shape[0]) +
@@ -2350,7 +2391,7 @@ def dotarg(__argname: str, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtyp
 
     Parameters
     ----------
-    argname : :class:`str`
+    name : :class:`str`
         The name of the argument.
     *arrays : :class:`Array` or something that can be :meth:`~Array.cast` into one
         The arrays to take inner products with.
@@ -2365,7 +2406,7 @@ def dotarg(__argname: str, *arrays: IntoArray, shape: Tuple[int, ...] = (), dtyp
         The inner product with shape ``shape + arrays[0].shape[1:] + ... + arrays[-1].shape[1:]``.
     '''
 
-    result = Argument(__argname, tuple(array.shape[0] for array in arrays) + tuple(shape), dtype=dtype)
+    result = Argument(name, tuple(array.shape[0] for array in arrays) + tuple(shape), dtype=dtype)
     for array in arrays:
         result = numpy.sum(_append_axes(result.transpose((*range(1, result.ndim), 0)), array.shape[1:]) * array, result.ndim-1)
     return result
@@ -2376,12 +2417,36 @@ class factor(Array):
 
     def __init__(self, array: Array) -> None:
         self._array = evaluable.factor(array)
-        super().__init__(shape=array.shape, dtype=array.dtype, spaces=set(), arguments=array.arguments)
+        super().__init__(shape=array.shape, dtype=array.dtype, spaces=set(), arguments=array._arguments)
 
     def lower(self, args: LowerArgs) -> evaluable.Array:
         return evaluable.prependaxes(self._array, args.points_shape)
 
 
+@nutils_dispatch
+def arguments_for(*arrays) -> Dict[str, Argument]:
+    '''Get all arguments that array(s) depend on.
+
+    Given any number of arrays, return a dictionary of all arguments involved,
+    mapping the name to the :class:`Argument` object. Raise a ``ValueError`` if
+    arrays have conflicting arguments, i.e. sharing a name but differing in
+    shape and/or dtype.
+    '''
+
+    arguments = {}
+    for array in arrays:
+        if isinstance(array, Array):
+            for name, (shape, dtype) in array._arguments.items():
+                argument = arguments.get(name)
+                if argument is None:
+                    arguments[name] = Argument(name, shape, dtype)
+                elif argument.shape != shape:
+                    raise ValueError(f'inconsistent shapes for argument {name!r}')
+                elif argument.dtype != dtype:
+                    raise ValueError(f'inconsistent dtypes for argument {name!r}')
+    return arguments
+
+
 # BASES
 
 
@@ -2443,7 +2508,7 @@ def __init__(self, ndofs: int, nelems: int, index: Array, coords: Array) -> None
         self.nelems = nelems
         self.index = Array.cast(index, dtype=int, ndim=0)
         self.coords = coords
-        arguments = _join_arguments((index.arguments, coords.arguments))
+        arguments = _join_arguments((index._arguments, coords._arguments))
         super().__init__((ndofs,), float, spaces=index.spaces | coords.spaces, arguments=arguments)
 
         _index = evaluable.Argument('_index', shape=(), dtype=int)
diff --git a/nutils/mesh.py b/nutils/mesh.py
index cb85a8e55..78470dac7 100644
--- a/nutils/mesh.py
+++ b/nutils/mesh.py
@@ -835,7 +835,7 @@ def unitcircle(nelems: int, variant: str) -> Tuple[Topology, function.Array]:
         # denominator onto a second order basis for efficient evaluation and
         # correct boundary gradients at the patch interfaces.
 
-        basis = topo.basis('spline', 2)
+        basis = topo.basis('spline', degree=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)
diff --git a/nutils/sample.py b/nutils/sample.py
index ca992e8af..bd76ed8fb 100644
--- a/nutils/sample.py
+++ b/nutils/sample.py
@@ -913,7 +913,7 @@ class _Integral(function.Array):
     def __init__(self, integrand: function.Array, sample: Sample) -> None:
         self._integrand = integrand
         self._sample = sample
-        super().__init__(shape=integrand.shape, dtype=float if integrand.dtype in (bool, int) else integrand.dtype, spaces=integrand.spaces - frozenset(sample.spaces), arguments=integrand.arguments)
+        super().__init__(shape=integrand.shape, dtype=float if integrand.dtype in (bool, int) else integrand.dtype, spaces=integrand.spaces - frozenset(sample.spaces), arguments=integrand._arguments)
 
     def lower(self, args: function.LowerArgs) -> evaluable.Array:
         ielem = evaluable.loop_index('_sample_' + '_'.join(self._sample.spaces), self._sample.nelems)
@@ -928,7 +928,7 @@ class _ConcatenatePoints(function.Array):
     def __init__(self, func: function.Array, sample: _TransformChainsSample) -> None:
         self._func = func
         self._sample = sample
-        super().__init__(shape=(sample.npoints, *func.shape), dtype=func.dtype, spaces=func.spaces - frozenset(sample.spaces), arguments=func.arguments)
+        super().__init__(shape=(sample.npoints, *func.shape), dtype=func.dtype, spaces=func.spaces - frozenset(sample.spaces), arguments=func._arguments)
 
     def lower(self, args: function.LowerArgs) -> evaluable.Array:
         axis = len(args.points_shape)
@@ -948,7 +948,7 @@ def __init__(self, func: function.Array, indices: evaluable.Array) -> None:
         self._func = func
         self._indices = indices
         assert indices.ndim == 1 and func.shape[0] == indices.shape[0].__index__()
-        super().__init__(shape=func.shape, dtype=func.dtype, spaces=func.spaces, arguments=func.arguments)
+        super().__init__(shape=func.shape, dtype=func.dtype, spaces=func.spaces, arguments=func._arguments)
 
     def lower(self, args: function.LowerArgs) -> evaluable.Array:
         func = self._func.lower(args)
diff --git a/nutils/solver.py b/nutils/solver.py
index beafcb65e..f4fb54e3e 100644
--- a/nutils/solver.py
+++ b/nutils/solver.py
@@ -1002,9 +1002,8 @@ def _split_trial_test(target):
 def _target_helper(target, *args):
     trial, test = _split_trial_test(target)
     if test is not None:
-        arguments = function._join_arguments(arg.arguments for arg in args)
-        testargs = [function.Argument(t, *arguments[t]) for t in test]
-        args = [map(arg.derivative, testargs) for arg in args]
+        arguments = function.arguments_for(*args)
+        args = [[arg.derivative(arguments[t]) for t in test] for arg in args]
     elif len(args) > 1:
         shapes = [{f.shape for f in ziparg if f is not None} for ziparg in zip(*args)]
         if any(len(arg) != len(shapes) for arg in args) or any(len(shape) != 1 for shape in shapes):
diff --git a/nutils/topology.py b/nutils/topology.py
index e6c7276c8..2776f06ab 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -341,19 +341,73 @@ def f_coords(self) -> function.Array:
 
         raise NotImplementedError
 
-    def basis(self, name: str, *args, **kwargs) -> function.Basis:
-        '''
-        Create a basis.
+    def _tensorial_bases(self, btype: str, **kwargs) -> function.Basis:
+        '''Generate bases for different tensorial directions as appropriate.'''
+
+        # This default implementation generates a single basis through dispatch
+        # to the basis_<btype> method, after stripping any hierarchical prefix
+        # to allow the same basis type to be used on the base topology as on
+        # its refinements.
+
+        if not self.ndims:
+            return
+
+        if btype.startswith('h-'):
+            btype = btype[2:]
+        elif btype.startswith('th-'):
+            btype = btype[3:]
+            kwargs.pop('truncation_tolerance', None)
+
+        f = getattr(self, 'basis_' + btype)
+        yield f(**kwargs)
+
+    @log.withcontext
+    def basis(self, btype=None, __degree=None, **kwargs):
+        '''Create a basis.
+
+        A basis is a set of basis functions organized in a one dimensional
+        array. The available types varies depending on the topology.
         '''
+
+        # The basis method is not redefined in Topology subclasses;
+        # specialization happens via the _tensorial_bases method.
+
+        if btype is None:
+            if 'name' not in kwargs:
+                raise TypeError("basis() missing 1 required positional argument: 'btype'")
+            btype = kwargs.pop('name')
+            warnings.deprecation("the basis type parameter has been renamed from 'name' to 'btype', please update your code")
+
+        if __degree is not None:
+            if 'degree' in kwargs:
+                raise TypeError("basis() got multiple values for argument 'degree'")
+            kwargs['degree'] = __degree
+            warnings.deprecation('the degree parameter must be specified as a keyword argument, please update your code')
+
         if self.ndims == 0:
             return function.PlainBasis([[[1]]], [[0]], 1, self.f_index, self.f_coords)
-        split = name.split('-', 1)
-        if len(split) == 2 and split[0] in ('h', 'th'):
-            name = split[1]  # default to non-hierarchical bases
-            if split[0] == 'th':
-                kwargs.pop('truncation_tolerance', None)
-        f = getattr(self, 'basis_' + name)
-        return f(*args, **kwargs)
+
+        basis, *other_bases = self._tensorial_bases(btype, **kwargs)
+        for other_basis in other_bases:
+            basis = numpy.ravel(basis[..., numpy.newaxis] * other_basis)
+        return basis
+
+    @log.withcontext
+    def field(self, name: str, /, *, shape=(), dtype=float, **kwargs):
+        '''Create a field.
+
+        A field is a contraction of a basis with an argument of given type and
+        shape, the latter prefixed with the required dof axes. The following
+        two functions are fully equivalent on non-tensorial topologies, and
+        functionally equivalent on all topologies:
+
+        ```
+        self.field(name, shape=S, dtype=D, **B)
+        function.field(name, self.basis(**B), shape=S, dtype=D)
+        ```
+        '''
+
+        return function.field(name, *self._tensorial_bases(**kwargs), shape=shape, dtype=dtype)
 
     def sample(self, ischeme: str, degree: int) -> Sample:
         'Create sample.'
@@ -1265,8 +1319,7 @@ def interfaces_spaces_unchecked(self, spaces: FrozenSet[str]) -> Topology:
             interfaces.append(self.topo1.interfaces_spaces(spaces1) * self.topo2)
         return Topology.disjoint_union(*interfaces)
 
-    def basis(self, name: str, degree: Union[int, Sequence[int]], **kwargs) -> function.Basis:
-        kwargs['degree'] = degree
+    def _tensorial_bases(self, btype, **kwargs):
         kwargs1 = {}
         kwargs2 = {}
 
@@ -1308,10 +1361,8 @@ def basis(self, name: str, degree: Union[int, Sequence[int]], **kwargs) -> funct
         kwargs1.update(kwargs)
         kwargs2.update(kwargs)
 
-        basis1 = self.topo1.basis(name, **kwargs1)
-        basis2 = self.topo2.basis(name, **kwargs2)
-        assert basis1.ndim == basis2.ndim == 1
-        return numpy.ravel(basis1[:,None] * basis2[None,:])
+        yield from self.topo1._tensorial_bases(btype, **kwargs1)
+        yield from self.topo2._tensorial_bases(btype, **kwargs2)
 
     def sample(self, ischeme: str, degree: int) -> Sample:
         return self.topo1.sample(ischeme, degree) * self.topo2.sample(ischeme, degree)
@@ -1389,8 +1440,8 @@ def f_coords(self) -> function.Array:
     def connectivity(self) -> Sequence[Sequence[int]]:
         return self.parent.connectivity
 
-    def basis(self, name: str, *args, **kwargs) -> function.Basis:
-        return self.parent.basis(name, *args, **kwargs)
+    def _tensorial_bases(self, btype, **kwargs) -> function.Basis:
+        yield from self.parent._tensorial_bases(btype, **kwargs)
 
     def sample(self, ischeme: str, degree: int) -> Sample:
         return self.parent.sample(ischeme, degree)
@@ -1851,8 +1902,8 @@ def points(self):
             topo = topo.basetopo
         return UnionTopology(ptopos, pnames)
 
-    def basis(self, name, *args, **kwargs):
-        return self.basetopo.basis(name, *args, **kwargs)
+    def _tensorial_bases(self, btype, **kwargs):
+        yield from self.basetopo._tensorial_bases(btype, **kwargs)
 
     @cached_property
     def refined(self):
@@ -2685,12 +2736,11 @@ def interfaces(self):
             irefs[iielem] = ref
         return SubsetTopology(baseinterfaces, irefs)
 
-    @log.withcontext
-    def basis(self, name, *args, **kwargs):
+    def _tensorial_bases(self, btype, **kwargs):
         if isinstance(self.basetopo, HierarchicalTopology):
             warnings.warn('basis may be linearly dependent; a linearly indepent basis is obtained by trimming first, then creating hierarchical refinements')
-        basis = self.basetopo.basis(name, *args, **kwargs)
-        return function.PrunedBasis(basis, self._indices, self.f_index, self.f_coords)
+        basis = self.basetopo.basis(btype, **kwargs)
+        yield function.PrunedBasis(basis, self._indices, self.f_index, self.f_coords)
 
     def locate(self, geom, coords, *, eps=0, skip_missing=False, **kwargs):
         sample = self.basetopo.locate(geom, coords, eps=eps, skip_missing=skip_missing, **kwargs)
@@ -2903,8 +2953,7 @@ def interfaces(self):
                 hopposites.append(level.interfaces.opposites[selection])
         return TransformChainsTopology(self.space, hreferences, transformseq.chain(htransforms, self.transforms.todims, self.ndims-1), transformseq.chain(hopposites, self.transforms.todims, self.ndims-1))
 
-    @log.withcontext
-    def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs):
+    def _tensorial_bases(self, btype, truncation_tolerance=1e-15, **kwargs):
         '''Create hierarchical basis.
 
         A hierarchical basis is constructed from bases on different levels of
@@ -2945,14 +2994,15 @@ def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs):
         basis : :class:`nutils.function.Array`
         '''
 
-        if name.startswith('h-'):
+        if btype.startswith('h-'):
             truncated = False
-            name = name[2:]
-        elif name.startswith('th-'):
+            btype = btype[2:]
+        elif btype.startswith('th-'):
             truncated = True
-            name = name[3:]
+            btype = btype[3:]
         else:
-            return super().basis(name, *args, **kwargs)
+            yield from super()._tensorial_bases(btype, **kwargs)
+            return
 
         # 1. identify active (supported) and passive (unsupported) basis functions
         ubases = []
@@ -2971,7 +3021,7 @@ def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs):
                 prev_ielems = ielems_i = numpy.unique(numpy.concatenate([numpy.asarray(touchielems_i, dtype=int), nontouchielems_i], axis=0))
                 prev_transforms = topo.transforms
 
-                basis_i = topo.basis(name, *args, **kwargs)
+                basis_i = topo.basis(btype, **kwargs)
                 assert isinstance(basis_i, function.Basis)
                 ubases.insert(0, basis_i)
                 # Basis functions that have at least one touchelem in their support.
@@ -3050,7 +3100,7 @@ def basis(self, name, *args, truncation_tolerance=1e-15, **kwargs):
                 degree = poly.degree(self.ndims, max(c.shape[-1] for c in trans_coeffs))
                 hbasis_coeffs.append(numpy.concatenate([poly.change_degree(c, self.ndims, degree) for c in trans_coeffs], axis=0))
 
-        return function.PlainBasis(hbasis_coeffs, hbasis_dofs, ndofs, self.f_index, self.f_coords)
+        yield function.PlainBasis(hbasis_coeffs, hbasis_dofs, ndofs, self.f_index, self.f_coords)
 
 
 class MultipatchTopology(TransformChainsTopology):
diff --git a/tests/test_evaluable.py b/tests/test_evaluable.py
index 338bb7784..7486e78a9 100644
--- a/tests/test_evaluable.py
+++ b/tests/test_evaluable.py
@@ -559,6 +559,7 @@ def _check(name, op, n_op, *arg_values, hasgrad=True, zerograd=False, ndim=2):
 _check('inflate-diagonal', lambda f: evaluable.Inflate(evaluable.Inflate(f, evaluable.constant(1), evaluable.constant(3)), evaluable.constant(1), evaluable.constant(3)), lambda a: numpy.diag(numpy.array([0, a, 0])), numpy.array(.5))
 _check('inflate-one', lambda f: evaluable.Inflate(f, evaluable.constant(0), evaluable.constant(1)), lambda a: numpy.array([a]), numpy.array(.5))
 _check('inflate-range', lambda f: evaluable.Inflate(f, evaluable.Range(evaluable.constant(3)), evaluable.constant(3)), lambda a: a, ANY(3))
+_check('inflate-twice', lambda f: evaluable._inflate(evaluable.Inflate(f, evaluable.constant([[0,1],[2,0]]), evaluable.constant(3)), evaluable.constant([0, 1, 0]), evaluable.constant(2), axis=0), lambda a: numpy.einsum('ij,jklm,lmn->ikn', [[1,0,1],[0,1,0]], a, [[[1,0,0],[0,1,0]],[[0,0,1],[1,0,0]]]), ANY(3,5,2,2))
 _check('take', lambda f: evaluable.Take(f, evaluable.constant([0, 3, 2])), lambda a: a[:, [0, 3, 2]], ANY(2, 4))
 _check('take-duplicate', lambda f: evaluable.Take(f, evaluable.constant([0, 3, 0])), lambda a: a[:, [0, 3, 0]], ANY(2, 4))
 _check('choose', lambda a, b, c: evaluable.Choose(a % 2, evaluable.stack([b, c], -1)), lambda a, b, c: numpy.choose(a % 2, [b, c]), INT(3, 3), ANY(3, 3), ANY(3, 3))
diff --git a/tests/test_expression_v2.py b/tests/test_expression_v2.py
index 437170439..91f77345e 100644
--- a/tests/test_expression_v2.py
+++ b/tests/test_expression_v2.py
@@ -560,27 +560,27 @@ def test_define_for_3d(self):
     def test_add_single_field(self):
         ns = expression_v2.Namespace()
         ns.add_field('u', numpy.array([1,2,3]))
-        self.assertEqual(ns.u.argshapes, dict(u=(3,)))
+        self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,)))
         self.assertEqual(ns.u.shape, ())
 
     def test_add_multiple_fields(self):
         ns = expression_v2.Namespace()
         ns.add_field(('u', 'v'), numpy.array([1,2,3]))
-        self.assertEqual(ns.u.argshapes, dict(u=(3,)))
+        self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,)))
         self.assertEqual(ns.u.shape, ())
-        self.assertEqual(ns.v.argshapes, dict(v=(3,)))
+        self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.v).values()}, dict(v=(3,)))
         self.assertEqual(ns.v.shape, ())
 
     def test_add_single_field_multiple_bases(self):
         ns = expression_v2.Namespace()
         ns.add_field('u', numpy.array([1,2,3]), numpy.array([4,5,6,7]))
-        self.assertEqual(ns.u.argshapes, dict(u=(3,4)))
+        self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,4)))
         self.assertEqual(ns.u.shape, ())
 
     def test_add_single_field_with_shape(self):
         ns = expression_v2.Namespace()
         ns.add_field('u', numpy.array([1,2,3]), shape=(2,))
-        self.assertEqual(ns.u.argshapes, dict(u=(3,2)))
+        self.assertEqual({a.name: a.shape for a in function.arguments_for(ns.u).values()}, dict(u=(3,2)))
         self.assertEqual(ns.u.shape, (2,))
 
     def test_copy(self):
diff --git a/tests/test_function.py b/tests/test_function.py
index 4fe6bb486..ba1bfb02c 100644
--- a/tests/test_function.py
+++ b/tests/test_function.py
@@ -128,7 +128,7 @@ def test_argshapes(self):
         a = function.Argument('a', (2, 3), dtype=int)
         b = function.Argument('b', (3,), dtype=int)
         f = (a * b[None]).sum(-1)
-        self.assertEqual(dict(f.argshapes), dict(a=(2, 3), b=(3,)))
+        self.assertEqual({a.name: a.shape for a in function.arguments_for(f).values()}, dict(a=(2, 3), b=(3,)))
 
     def test_argshapes_shape_mismatch(self):
         with self.assertRaises(Exception):
@@ -425,7 +425,7 @@ def test(self):
         f = function._Unlower(e, frozenset(), arguments, function.LowerArgs((2, 3), {}, {}))
         self.assertEqual(f.shape, (4, 5))
         self.assertEqual(f.dtype, int)
-        self.assertEqual(f.arguments, arguments)
+        self.assertEqual(f._arguments, arguments)
         self.assertEqual(f.lower(function.LowerArgs((2, 3), {}, {})), e)
         with self.assertRaises(ValueError):
             f.lower(function.LowerArgs((3, 4), {}, {}))
@@ -790,11 +790,11 @@ def test_different_dtype(self):
 
     def test_nonempty_spaces(self):
         topo, geom = mesh.unitsquare(1, 'square')
-        with self.assertRaisesRegex(ValueError, "replacement functions cannot contain spaces, but replacement for Argument 'foo' contains space X."):
+        with self.assertRaisesRegex(ValueError, "replacement functions cannot be bound to a space, but replacement for Argument 'foo' is bound to X."):
             function.replace_arguments(function.Argument('foo', (2,), dtype=float), dict(foo=geom))
 
 
-class dotarg(TestCase):
+class field(TestCase):
 
     def assertEvalAlmostEqual(self, f_actual, desired, **arguments):
         self.assertEqual(f_actual.shape, desired.shape)
@@ -810,13 +810,13 @@ def test(self):
         a243 = numpy.arange(24, dtype=float).reshape(2, 4, 3)
         a4 = numpy.arange(4, dtype=float)
         a45 = numpy.arange(20, dtype=float).reshape(4, 5)
-        self.assertEvalAlmostEqual(function.dotarg('arg'), a, arg=a)
-        self.assertEvalAlmostEqual(function.dotarg('arg', shape=(2, 3)), a23, arg=a23)
-        self.assertEvalAlmostEqual(function.dotarg('arg', a4), numpy.einsum('i,i->', a4, a4), arg=a4)
-        self.assertEvalAlmostEqual(function.dotarg('arg', a45), numpy.einsum('i,ij->j', a4, a45), arg=a4)
-        self.assertEvalAlmostEqual(function.dotarg('arg', a24, shape=(3,)), numpy.einsum('ij,ik->jk', a23, a24), arg=a23)
-        self.assertEvalAlmostEqual(function.dotarg('arg', a2, a45), numpy.einsum('ij,i,jk->k', a24, a2, a45), arg=a24)
-        self.assertEvalAlmostEqual(function.dotarg('arg', a2, a45, shape=(3,)), numpy.einsum('ijk,i,jl->kl', a243, a2, a45), arg=a243)
+        self.assertEvalAlmostEqual(function.field('arg'), a, arg=a)
+        self.assertEvalAlmostEqual(function.field('arg', shape=(2, 3)), a23, arg=a23)
+        self.assertEvalAlmostEqual(function.field('arg', a4), numpy.einsum('i,i->', a4, a4), arg=a4)
+        self.assertEvalAlmostEqual(function.field('arg', a45), numpy.einsum('i,ij->j', a4, a45), arg=a4)
+        self.assertEvalAlmostEqual(function.field('arg', a24, shape=(3,)), numpy.einsum('ij,ik->jk', a23, a24), arg=a23)
+        self.assertEvalAlmostEqual(function.field('arg', a2, a45), numpy.einsum('ij,i,jk->k', a24, a2, a45), arg=a24)
+        self.assertEvalAlmostEqual(function.field('arg', a2, a45, shape=(3,)), numpy.einsum('ijk,i,jl->kl', a243, a2, a45), arg=a243)
 
 
 class jacobian(TestCase):
@@ -1388,13 +1388,13 @@ def test_stokes(self):
 class Eval(TestCase):
 
     def test_single(self):
-        f = function.dotarg('v', numpy.array([1, 2, 3]))
+        f = function.field('v', numpy.array([1, 2, 3]))
         retval = function.eval(f, v=numpy.array([4, 5, 6]))
         self.assertEqual(retval, 4+10+18)
 
     def test_multiple(self):
-        f = function.dotarg('v', numpy.array([1, 2, 3]))
-        g = function.dotarg('v', numpy.array([3, 2, 1]))
+        f = function.field('v', numpy.array([1, 2, 3]))
+        g = function.field('v', numpy.array([3, 2, 1]))
         retvals = function.eval([f, g], v=numpy.array([4, 5, 6]))
         self.assertEqual(retvals, (4+10+18, 12+10+6))
 
@@ -1433,13 +1433,13 @@ class factor(TestCase):
 
     def test_lower_with_points(self):
         topo, geom = mesh.rectilinear([3])
-        f = function.dotarg('dof')
+        f = function.field('dof')
         v = topo.sample('uniform', 1).eval(function.factor(f**2) * geom[0], dof=2.)
         self.assertAllAlmostEqual(v, [2, 6, 10])
 
     def test_constant(self):
         topo, geom = mesh.rectilinear([3])
-        basis = topo.basis('std', 1)
+        basis = topo.basis('std', degree=1)
         f = function.factor(topo.integral(basis[:, numpy.newaxis] * basis, degree=2))
         ((i, j), v, shape), = function.evaluate(f, _post=sparse.extract)
         self.assertAllEqual(i, [0, 0, 1, 1, 1, 2, 2, 2, 3, 3])
@@ -1450,3 +1450,30 @@ def test_lower_spaces(self):
         topo, geom = mesh.rectilinear([3])
         with self.assertRaisesRegex(ValueError, 'cannot lower function with spaces \(.+\) - did you forget integral or sample?'):
             function.factor(geom)
+
+
+class arguments_for(TestCase):
+
+    def test_single(self):
+        x = function.field('x', numpy.array([1,2,3]), shape=(2,), dtype=int)
+        f = x**2
+        self.assertEqual({a.name: (a.shape, a.dtype) for a in function.arguments_for(f).values()},
+            {'x': ((3,2), int)})
+
+    def test_multiple(self):
+        x = function.field('x', numpy.array([1,2,3]), shape=(2,), dtype=int)
+        y = function.field('y', numpy.array([4,5]), dtype=float)
+        z = function.field('z', dtype=complex)
+        f = x * y
+        g = x**2 * z
+        self.assertEqual({a.name: (a.shape, a.dtype) for a in function.arguments_for(f, g).values()},
+            {'x': ((3,2), int), 'y': ((2,), float), 'z': ((), complex)})
+
+    def test_conflict(self):
+        x = function.field('x', numpy.array([1,2,3]), shape=(2,), dtype=int)
+        y1 = function.field('y', numpy.array([4,5]), dtype=float)
+        y2 = function.field('y', dtype=complex)
+        f = x * y1
+        g = x**2 * y2
+        with self.assertRaisesRegex(ValueError, "inconsistent shapes for argument 'y'"):
+            function.arguments_for(f, g)
diff --git a/tests/test_solver.py b/tests/test_solver.py
index 13f7b5822..546c5ac00 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -101,8 +101,8 @@ def setUp(self):
         super().setUp()
         domain, geom = mesh.rectilinear([8, 8])
         basis = domain.basis('std', degree=1)
-        u = function.dotarg('u', basis)
-        v = function.dotarg('v', basis)
+        u = function.field('u', basis)
+        v = function.field('v', basis)
         sqr = domain.boundary['left'].integral(u**2, degree=2)
         self.cons = solver.optimize('u,', sqr)
         self.residual = domain.integral((v.grad(geom) @ u.grad(geom))*function.J(geom), degree=2) \
@@ -205,9 +205,9 @@ def setUp(self):
         domain, geom = mesh.rectilinear([numpy.linspace(0, 1, 9)] * 2)
         ubasis = domain.basis('std', degree=2)
         if self.vector:
-            u = function.dotarg('dofs', ubasis.vector(2))
+            u = function.field('dofs', ubasis.vector(2))
         else:
-            u = function.dotarg('dofs', ubasis, shape=(2,))
+            u = function.field('dofs', ubasis, shape=(2,))
         Geom = geom * [1.1, 1] + u
         self.cons = solver.optimize('dofs', domain.boundary['left,right'].integral((u**2).sum(0), degree=4), droptol=1e-15)
         self.boolcons = ~numpy.isnan(self.cons)
@@ -258,7 +258,7 @@ def setUp(self):
         self.domain, geom = mesh.rectilinear([2, 2])
         self.ubasis = self.domain.basis('std', degree=1)
         self.ns = Namespace()
-        self.ns.u = function.dotarg('dofs', self.ubasis)
+        self.ns.u = function.field('dofs', self.ubasis)
 
     def test_linear(self):
         err = self.domain.boundary['bottom'].integral('(u - 1)^2' @ self.ns, degree=2)
@@ -337,8 +337,8 @@ class System(TestCase):
 
     def test_constant(self):
         domain, geom = mesh.rectilinear([9])
-        u = function.dotarg('u', domain.basis('std', 1))
-        v = function.dotarg('v', domain.basis('std', 1))
+        u = function.field('u', domain.basis('std', degree=1))
+        v = function.field('v', domain.basis('std', degree=1))
         f = domain.integral(v * (1 + u) * function.J(geom), degree=2)
         sys = solver.System(f, trial='u', test='v')
         self.assertTrue(sys.is_linear)
@@ -358,7 +358,7 @@ def test_constant(self):
 
     def test_constant_symmetric(self):
         domain, geom = mesh.rectilinear([9])
-        u = function.dotarg('u', domain.basis('std', 1))
+        u = function.field('u', domain.basis('std', degree=1))
         f = domain.integral((1 + u - u**2) * function.J(geom), degree=2)
         sys = solver.System(f, trial='u')
         self.assertTrue(sys.is_linear)
@@ -379,9 +379,9 @@ def test_constant_symmetric(self):
 
     def test_constant_matrix(self):
         domain, geom = mesh.rectilinear([9])
-        u = function.dotarg('u', domain.basis('std', 1))
-        v = function.dotarg('v', domain.basis('std', 1))
-        t = function.dotarg('t')
+        u = function.field('u', domain.basis('std', degree=1))
+        v = function.field('v', domain.basis('std', degree=1))
+        t = function.field('t')
         f = domain.integral(v * (t + u) * function.J(geom), degree=2)
         sys = solver.System(f, trial='u', test='v')
         self.assertTrue(sys.is_linear)
@@ -401,9 +401,9 @@ def test_constant_matrix(self):
 
     def test_linear(self):
         domain, geom = mesh.rectilinear([9])
-        u = function.dotarg('u', domain.basis('std', 1))
-        v = function.dotarg('v', domain.basis('std', 1))
-        t = function.dotarg('t')
+        u = function.field('u', domain.basis('std', degree=1))
+        v = function.field('v', domain.basis('std', degree=1))
+        t = function.field('t')
         f = domain.integral(v * (u * t + 1) * function.J(geom), degree=2)
         sys = solver.System(f, trial='u', test='v')
         self.assertTrue(sys.is_linear)
@@ -419,9 +419,9 @@ def test_linear(self):
 
     def test_nonlinear(self):
         domain, geom = mesh.rectilinear([9])
-        u = function.dotarg('u', domain.basis('std', 1))
-        v = function.dotarg('v', domain.basis('std', 1))
-        t = function.dotarg('t')
+        u = function.field('u', domain.basis('std', degree=1))
+        v = function.field('v', domain.basis('std', degree=1))
+        t = function.field('t')
         f = domain.integral(v * (u**2 + t) * function.J(geom), degree=2)
         sys = solver.System(f, trial='u', test='v')
         self.assertFalse(sys.is_linear)
@@ -437,7 +437,7 @@ def test_nonlinear(self):
 
     def test_nonlinear_symmetric(self):
         domain, geom = mesh.rectilinear([9])
-        u = function.dotarg('u', domain.basis('std', 1))
+        u = function.field('u', domain.basis('std', degree=1))
         f = domain.integral(numpy.exp(u) * function.J(geom), degree=2)
         sys = solver.System(f, trial='u')
         self.assertFalse(sys.is_linear)
@@ -457,7 +457,7 @@ class system_finitestrain(TestCase):
     def setUp(self):
         super().setUp()
         domain, geom = mesh.rectilinear([numpy.linspace(0, 1, 9)] * 2)
-        u = function.dotarg('u', domain.basis('std', degree=2), shape=(2,))
+        u = function.field('u', domain.basis('std', degree=2), shape=(2,))
         Geom = geom * [1.1, 1] + u
         self.cons = solver.optimize('u,', domain.boundary['left,right'].integral((u**2).sum(0), degree=4), droptol=1e-15)
         strain = .5 * (function.outer(Geom.grad(geom), axis=1).sum(0) - function.eye(2))
@@ -497,10 +497,10 @@ def setUp(self):
         dx = function.J(geom)
         ubasis = domain.basis('std', degree=2)
         pbasis = domain.basis('std', degree=1)
-        u = function.dotarg('u', ubasis, shape=(2,))
-        v = function.dotarg('v', ubasis, shape=(2,))
-        p = function.dotarg('p', pbasis)
-        q = function.dotarg('q', pbasis)
+        u = function.field('u', ubasis, shape=(2,))
+        v = function.field('v', ubasis, shape=(2,))
+        p = function.field('p', pbasis)
+        q = function.field('q', pbasis)
         self.cons = solver.optimize('u,', domain.boundary['top,bottom'].integral((u**2).sum(), degree=4), droptol=1e-10)
         self.cons = solver.optimize('u,', domain.boundary['left'].integral((u[0]-uin)**2 + u[1]**2, degree=4), droptol=1e-10, constrain=self.cons)
         res = gauss.integral((viscosity * numpy.einsum('ij,ij', v.grad(geom), u.grad(geom) + u.grad(geom).T) - v.div(geom) * p + q * u.div(geom)) * dx)
diff --git a/tests/test_topology.py b/tests/test_topology.py
index c1c9f428a..fb98db908 100644
--- a/tests/test_topology.py
+++ b/tests/test_topology.py
@@ -1190,7 +1190,7 @@ class multipatch_misc(TestCase):
     def test_basis_merge_dofs(self):
         patches = (0, 4, 2, 5), (9, 1, 5, 3), (2, 5, 6, 7), (5, 3, 7, 8)
         domain, geom = mesh.multipatch(patches=patches, nelems=1)
-        basis = domain.basis('spline', 1)
+        basis = domain.basis('spline', degree=1)
         actual = domain.sample('bezier', 2).eval(basis)
         desired = numpy.zeros((16, 10))
         desired[numpy.arange(16), [0, 1, 2, 3, 4, 5, 3, 6, 2, 3, 7, 8, 3, 6, 8, 9]] = 1
@@ -1207,7 +1207,7 @@ def test_ring(self):
         domain, geom = mesh.multipatch(patches=patches, patchverts=patchverts, nelems=1)
         self.assertEqual(domain.connectivity.tolist(), [[1, -1,  2, -1], [-1, 0, 2, -1], [1, -1, -1, 0]])
         self.assertEqual(len(domain.interfaces['interpatch']), 3)
-        self.assertEqual(len(domain.basis('std', 1)), 7)
+        self.assertEqual(len(domain.basis('std', degree=1)), 7)
         self.assertAlmostEqual(domain.integrate(function.J(geom), degree=1), 1)
 
     def test_partial_ring(self):
@@ -1221,7 +1221,7 @@ def test_partial_ring(self):
         domain, geom = mesh.multipatch(patches=patches, nelems=1)
         self.assertEqual(len(domain.interfaces['interpatch']), 5) # 5th interface via 15-12-5-1
         self.assertEqual(len(domain.boundary), 14)
-        self.assertEqual(len(domain.basis('std', 1)), 16)
+        self.assertEqual(len(domain.basis('std', degree=1)), 16)
         # this test fails if MultipatchTopology.basis_std delegates to
         # TransformChainsTopology.basis_std rather than MultipatchTopology.basis_spline.