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..5f218de5b 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'
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..d1d28f697 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(res.argshapes['p'], 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(sqr.argshapes['ψ'], 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 cda9bc8bf..ed86b0b32 100644
--- a/examples/platewithhole.py
+++ b/examples/platewithhole.py
@@ -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 be294675c..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.field('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/topology.py b/nutils/topology.py
index abbb2b93c..2776f06ab 100644
--- a/nutils/topology.py
+++ b/nutils/topology.py
@@ -392,6 +392,23 @@ def basis(self, btype=None, __degree=None, **kwargs):
             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.'