partition

docs/notes/binary_operations.rst

@@ -25,4 +25,4 @@ Binary Operations on Arrays
 
 1. In the above table the summation axes are numbered backward. For example, ``sum(-1)`` is used to sum over the last axis of an array. Although forward numbering is possible in many situations, backward numbering is generally preferred in Nutils code.
 2. When a summation over multiple axes is performed (#6), these axes are to be listed. In the case of single-axis summations listing is optional (for example ``sum(-1)`` is equivalent to ``sum([-1])``). The shorter notation ``sum(-1)`` is preferred.
-3. When the numer of dimensions of the two arguments of a binary operation mismatch, singleton axes are automatically prepended to the "shorter" argument. This property can be used to shorten notation. For example, #3 can be written as ``(A*b).sum(-1)``. To avoid ambiguities, in general, such abbreviations are discouraged.
+3. When the number of dimensions of the two arguments of a binary operation mismatch, singleton axes are automatically prepended to the "shorter" argument. This property can be used to shorten notation. For example, #3 can be written as ``(A*b).sum(-1)``. To avoid ambiguities, in general, such abbreviations are discouraged.

docs/tutorial.rst

@@ -47,8 +47,8 @@ Note that discretization inevitably implies approximation, i.e. :math:`u ≠
 space of piecewise linears, which contains the exact solution. We therefore
 expect our Finite Element solution to be exact.
 
-Wetting your appetite
----------------------
+Whetting your appetite
+----------------------
 
 The computation can be set up in about 20 lines of Nutils code, including
 visualization. The entire script is presented below, in copy-pasteable form
@@ -123,7 +123,7 @@ size between 0 and 1 is generated by
 
 .. console::
     >>> nutils.mesh.rectilinear([[0, 0.25, 0.5, 0.75, 1.0]])
-    (StructuredTopology<4>, Array<1>)
+    (StructuredLine<4>, Array<1>)
 
 Alternatively we could have used :func:`numpy.linspace` to generate a sequence
 of equidistant vertices, and unpack the resulting tuple:
@@ -142,7 +142,7 @@ is generated by
 
 .. console::
     >>> nutils.mesh.rectilinear([numpy.linspace(0, 1, 5), numpy.linspace(0, 1, 9)])
-    (StructuredTopology<4x8>, Array<2>)
+    (StructuredLine<4>*StructuredLine<8>, Array<2>)
 
 Any topology defines a boundary via the :attr:`Topology.boundary
 <nutils.topology.Topology.boundary>` attribute. Optionally, a topology can
@@ -152,7 +152,7 @@ dimension, making the left boundary accessible as:
 
 .. console::
     >>> topo.boundary['left']
-    StructuredTopology<>
+    PointsTopology<1>
 
 Optionally, a topology can be made periodic in one or more dimensions by
 passing a list of dimension indices to be periodic via the keyword argument
@@ -161,7 +161,7 @@ two-dimensional mesh periodic, add ``periodic=[1]``:
 
 .. console::
     >>> nutils.mesh.rectilinear([numpy.linspace(0, 1, 5), numpy.linspace(0, 1, 9)], periodic=[1])
-    (StructuredTopology<4x8p>, Array<2>)
+    (StructuredLine<4>*StructuredLine<8p>, Array<2>)
 
 Note that in this case the boundary topology, though still available, is empty.
 
@@ -177,7 +177,7 @@ it helps to think of a basis as evaluating always to the full array.
 
 Several :class:`~nutils.topology.Topology` objects support creating bases via
 the :meth:`Topology.basis() <nutils.topology.Topology.basis>` method.  A
-:class:`~nutils.topology.StructuredTopology`, as generated by
+:class:`~nutils.topology.StructuredLine`, as generated by
 :func:`nutils.mesh.rectilinear`, can create a spline basis with arbitrary
 degree and arbitrary continuity. The following generates a degree one spline
 basis on our previously created unit line topology ``topo``:
@@ -609,7 +609,7 @@ The optimization problem can also be solved by the
 .. console::
     >>> nutils.solver.optimize('lhs', sqr)
     optimize > solve > solving 5 dof system to machine precision using direct solver
-    optimize > solve > solver returned with residual 0e+00
+    optimize > solve > solver returned with residual 0e+00±1e-15
     optimize > optimum value 0.00e+00±1e-15
     array([0.  , 0.25, 0.5 , 0.75, 1.  ])±1e-15
 
@@ -676,7 +676,7 @@ second argument of :class:`Topology.sample()
 .. console::
     >>> bezier = topo.sample('bezier', 2)
     >>> bezier
-    Sample<1D, 4 elems, 8 points>
+    UniformSample<1D, 4 elems, 8 points>
 
 The resulting :class:`nutils.sample.Sample` object can be used to evaluate
 :class:`~nutils.function.Array` functions via the :meth:`Sample.eval(func)

examples/adaptivity.py

@@ -45,8 +45,8 @@ def main(etype:str, btype:str, degree:int, nrefine:int):
       if irefine:
         refdom = domain.refined
         ns.refbasis = refdom.basis(btype, degree=degree)
-        indicator = refdom.integral('refbasis_n,k u_,k d:x' @ ns, degree=degree*2).eval(lhs=lhs)
-        indicator -= refdom.boundary.integral('refbasis_n u_,k n_k d:x' @ ns, degree=degree*2).eval(lhs=lhs)
+        indicator = refdom.integral('d(refbasis_n, x_k) d(u, x_k) d:x' @ ns, degree=degree*2).eval(lhs=lhs)
+        indicator -= refdom.boundary.integral('refbasis_n d(u, x_k) n(x_k) d:x' @ ns, degree=degree*2).eval(lhs=lhs)
         supp = ns.refbasis.get_support(indicator**2 > numpy.mean(indicator**2))
         domain = domain.refined_by(ns.refbasis.transforms[supp])
 
@@ -62,11 +62,11 @@ def main(etype:str, btype:str, degree:int, nrefine:int):
       sqr = domain.boundary.integral('du^2 d:x' @ ns, degree=7)
       cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
 
-      res = domain.integral('basis_n,k u_,k d:x' @ ns, degree=degree*2)
+      res = domain.integral('d(basis_n, x_k) d(u, x_k) d:x' @ ns, degree=degree*2)
       lhs = solver.solve_linear('lhs', res, constrain=cons)
 
       ndofs = len(ns.basis)
-      error = domain.integral('<du^2, du_,k du_,k>_i d:x' @ ns, degree=7).eval(lhs=lhs)**.5
+      error = domain.integral('<du^2, d(du, x_k) d(du, x_k)>_i d:x' @ ns, degree=7).eval(lhs=lhs)**.5
       rate, offset = linreg.add(numpy.log(len(ns.basis)), numpy.log(error))
       treelog.user('ndofs: {ndofs}, L2 error: {error[0]:.2e} ({rate[0]:.2f}), H1 error: {error[1]:.2e} ({rate[1]:.2f})'.format(ndofs=len(ns.basis), error=error, rate=rate))
 

examples/burgers.py

@@ -36,8 +36,8 @@ def main(nelems:int, ndims:int, degree:int, timescale:float, newtontol:float, en
   ns.f = '.5 u^2'
   ns.C = 1
 
-  res = domain.integral('-basis_n,0 f d:x' @ ns, degree=5)
-  res += domain.interfaces.integral('-[basis_n] n_0 ({f} - .5 C [u] n_0) d:x' @ ns, degree=degree*2)
+  res = domain.integral('-d(basis_n, x_i) δ_i0 f d:x' @ ns, degree=5)
+  res += domain.interfaces.integral('-[basis_n] n(x_i) δ_i0 ({f} - .5 C [u] n(x_j) δ_j0) d:x' @ ns, degree=degree*2)
   inertia = domain.integral('basis_n u d:x' @ ns, degree=5)
 
   sqr = domain.integral('(u - exp(-?y_i ?y_i)(y_i = 5 (x_i - 0.5_i)))^2 d:x' @ ns, degree=5)
@@ -93,6 +93,6 @@ def test_1d_p2(self):
   def test_2d_p1(self):
     lhs = main(ndims=2, nelems=4, timescale=.1, degree=1, endtime=.01, newtontol=1e-5)
     self.assertAlmostEqual64(lhs, '''
-      eNoNyKENhEAQRuGEQsCv2SEzyQZHDbRACdsDJNsBjqBxSBxBHIgJ9xsqQJ1Drro1L1/eYBZceGz8njrR
-      yacm8UQLBvPYCw1airpyUVYSJLhKijK4IC01WDnqqxvX8OTl427aU73sctPGr3qqceBnRzOjo0xy9JpJ
-      R73m6R6YMZo/Q+FCLQ==''')
+      eNoNx6ENhjAQBtCEQcDX9Mhd0uCYgRUYoTtAwga4P2gcEkcQUHGBzzABCoes+slTrzcT4hk0aDwn9ObA
+      bQcKHHig2x6oUFOWF1JIltdUIerMnXTuIzNHfXVhL5vbnJeFXy3h6aJVVrnIU4kdj20okUQaeuyOnxmR
+      otVWU4zf/jlhQi0=''')

examples/cahnhilliard.py

@@ -72,9 +72,9 @@ def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[f
   ns.dt = timestep
 
   nrg_mix = domain.integral('F d:x' @ ns, degree=7)
-  nrg_iface = domain.integral('.5 c_,k c_,k d:x' @ ns, degree=7)
+  nrg_iface = domain.integral('.5 d(c, x_k) d(c, x_k) d:x' @ ns, degree=7)
   nrg_wall = domain.boundary.integral('(abs(ewall) + c ewall) d:x' @ ns, degree=7)
-  nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(dF - m dc - .5 dt epsilon^2 m_,k m_,k) d:x' @ ns, degree=7)
+  nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(dF - m dc - .5 dt epsilon^2 d(m, x_k) d(m, x_k)) d:x' @ ns, degree=7)
 
   numpy.random.seed(seed)
   lhs = numpy.random.normal(0, .5, ns.cbasis.shape) # initial condition

examples/cylinderflow.py

@@ -59,9 +59,9 @@ def main(nelems:int, degree:int, reynolds:float, rotation:float, timestep:float,
   ns.ubasis_ni = 'unbasis_n J_i0 + utbasis_n J_i1' # piola transformation
   ns.u_i = 'ubasis_ni ?lhs_n'
   ns.p = 'pbasis_n ?lhs_n'
-  ns.sigma_ij = '(u_i,j + u_j,i) / Re - p δ_ij'
+  ns.sigma_ij = '(d(u_i, x_j) + d(u_j, x_i)) / Re - p δ_ij'
   ns.N = 10 * degree / elemangle # Nitsche constant based on element size = elemangle/2
-  ns.nitsche_ni = '(N ubasis_ni - (ubasis_ni,j + ubasis_nj,i) n_j) / Re'
+  ns.nitsche_ni = '(N ubasis_ni - (d(ubasis_ni, x_j) + d(ubasis_nj, x_i)) n_j) / Re'
   ns.rotation = rotation
   ns.uwall_i = '0.5 rotation <-sin(phi), cos(phi)>_i'
 
@@ -75,7 +75,7 @@ def main(nelems:int, degree:int, reynolds:float, rotation:float, timestep:float,
   numpy.random.seed(seed)
   lhs0 *= numpy.random.normal(1, .1, lhs0.shape) # add small velocity noise
 
-  res = domain.integral('(ubasis_ni u_i,j u_j + ubasis_ni,j sigma_ij + pbasis_n u_k,k) d:x' @ ns, degree=9)
+  res = domain.integral('(ubasis_ni d(u_i, x_j) u_j + d(ubasis_ni, x_j) sigma_ij + pbasis_n d(u_k, x_k)) d:x' @ ns, degree=9)
   res += domain.boundary['inner'].integral('(nitsche_ni (u_i - uwall_i) - ubasis_ni sigma_ij n_j) d:x' @ ns, degree=9)
   inertia = domain.integral('ubasis_ni u_i d:x' @ ns, degree=9)
 

examples/drivencavity-compatible.py

@@ -42,18 +42,18 @@ def main(nelems:int, degree:int, reynolds:float):
   ns.u_i = 'ubasis_ni ?lhs_n'
   ns.p = 'pbasis_n ?lhs_n'
   ns.l = 'lbasis_n ?lhs_n'
-  ns.stress_ij = '(u_i,j + u_j,i) / Re - p δ_ij'
+  ns.stress_ij = '(d(u_i, x_j) + d(u_j, x_i)) / Re - p δ_ij'
   ns.uwall = domain.boundary.indicator('top'), 0
   ns.N = 5 * degree * nelems # nitsche constant based on element size = 1/nelems
-  ns.nitsche_ni = '(N ubasis_ni - (ubasis_ni,j + ubasis_nj,i) n_j) / Re'
+  ns.nitsche_ni = '(N ubasis_ni - (d(ubasis_ni, x_j) + d(ubasis_nj, x_i)) n(x_j)) / Re'
 
-  res = domain.integral('(ubasis_ni,j stress_ij + pbasis_n (u_k,k + l) + lbasis_n p) d:x' @ ns, degree=2*degree)
+  res = domain.integral('(d(ubasis_ni, x_j) stress_ij + pbasis_n (d(u_k, x_k) + l) + lbasis_n p) d:x' @ ns, degree=2*degree)
   res += domain.boundary.integral('(nitsche_ni (u_i - uwall_i) - ubasis_ni stress_ij n_j) d:x' @ ns, degree=2*degree)
   with treelog.context('stokes'):
     lhs0 = solver.solve_linear('lhs', res)
     postprocess(domain, ns, lhs=lhs0)
 
-  res += domain.integral('ubasis_ni u_i,j u_j d:x' @ ns, degree=3*degree)
+  res += domain.integral('ubasis_ni d(u_i, x_j) u_j d:x' @ ns, degree=3*degree)
   with treelog.context('navierstokes'):
     lhs1 = solver.newton('lhs', res, lhs0=lhs0).solve(tol=1e-10)
     postprocess(domain, ns, lhs=lhs1)
@@ -66,13 +66,14 @@ def main(nelems:int, degree:int, reynolds:float):
 
 def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
 
-  div = domain.integral('(u_k,k)^2 d:x' @ ns, degree=1).eval(**arguments)**.5
+  div = domain.integral('d(u_k, x_k)^2 d:x' @ ns, degree=1).eval(**arguments)**.5
   treelog.info('velocity divergence: {:.2e}'.format(div)) # confirm that velocity is pointwise divergence-free
 
   ns = ns.copy_() # copy namespace so that we don't modify the calling argument
   ns.streambasis = domain.basis('std', degree=2)[1:] # remove first dof to obtain non-singular system
   ns.stream = 'streambasis_n ?streamdofs_n' # stream function
-  sqr = domain.integral('((u_0 - stream_,1)^2 + (u_1 + stream_,0)^2) d:x' @ ns, degree=4)
+  ns.ε = function.levicivita(2)
+  sqr = domain.integral('(u_i - ε_ij d(stream, x_j)) (u_i - ε_ij d(stream, x_j)) d:x' @ ns, degree=4)
   arguments['streamdofs'] = solver.optimize('streamdofs', sqr, arguments=arguments) # compute streamlines
 
   bezier = domain.sample('bezier', 9)

examples/drivencavity.py

@@ -35,7 +35,7 @@ def main(nelems:int, etype:str, degree:int, reynolds:float):
   ])
   ns.u_i = 'ubasis_ni ?lhs_n'
   ns.p = 'pbasis_n ?lhs_n'
-  ns.stress_ij = '(u_i,j + u_j,i) / Re - p δ_ij'
+  ns.stress_ij = '(d(u_i, x_j) + d(u_j, x_i)) / Re - p δ_ij'
 
   sqr = domain.boundary.integral('u_k u_k d:x' @ ns, degree=degree*2)
   wallcons = solver.optimize('lhs', sqr, droptol=1e-15)
@@ -46,12 +46,12 @@ def main(nelems:int, etype:str, degree:int, reynolds:float):
   cons = numpy.choose(numpy.isnan(lidcons), [lidcons, wallcons])
   cons[-1] = 0 # pressure point constraint
 
-  res = domain.integral('(ubasis_ni,j stress_ij + pbasis_n u_k,k) d:x' @ ns, degree=degree*2)
+  res = domain.integral('(d(ubasis_ni, x_j) stress_ij + pbasis_n d(u_k, x_k)) d:x' @ ns, degree=degree*2)
   with treelog.context('stokes'):
     lhs0 = solver.solve_linear('lhs', res, constrain=cons)
     postprocess(domain, ns, lhs=lhs0)
 
-  res += domain.integral('.5 (ubasis_ni u_i,j - ubasis_ni,j u_i) u_j d:x' @ ns, degree=degree*3)
+  res += domain.integral('.5 (ubasis_ni d(u_i, x_j) - d(ubasis_ni, x_j) u_i) u_j d:x' @ ns, degree=degree*3)
   with treelog.context('navierstokes'):
     lhs1 = solver.newton('lhs', res, lhs0=lhs0, constrain=cons).solve(tol=1e-10)
     postprocess(domain, ns, lhs=lhs1)
@@ -67,7 +67,8 @@ def postprocess(domain, ns, every=.05, spacing=.01, **arguments):
   ns = ns.copy_() # copy namespace so that we don't modify the calling argument
   ns.streambasis = domain.basis('std', degree=2)[1:] # remove first dof to obtain non-singular system
   ns.stream = 'streambasis_n ?streamdofs_n' # stream function
-  sqr = domain.integral('((u_0 - stream_,1)^2 + (u_1 + stream_,0)^2) d:x' @ ns, degree=4)
+  ns.ε = function.levicivita(2)
+  sqr = domain.integral('(u_i - ε_ij d(stream, x_j)) (u_i - ε_ij d(stream, x_j)) d:x' @ ns, degree=4)
   arguments['streamdofs'] = solver.optimize('streamdofs', sqr, arguments=arguments) # compute streamlines
 
   bezier = domain.sample('bezier', 9)

examples/elasticity.py

@@ -34,14 +34,14 @@ def main(nelems:int, etype:str, btype:str, degree:int, poisson:float):
   ns.X_i = 'x_i + u_i'
   ns.lmbda = 2 * poisson
   ns.mu = 1 - 2 * poisson
-  ns.strain_ij = '(u_i,j + u_j,i) / 2'
+  ns.strain_ij = '(d(u_i, x_j) + d(u_j, x_i)) / 2'
   ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
 
   sqr = domain.boundary['left'].integral('u_k u_k d:x' @ ns, degree=degree*2)
   sqr += domain.boundary['right'].integral('(u_0 - .5)^2 d:x' @ ns, degree=degree*2)
   cons = solver.optimize('lhs', sqr, droptol=1e-15)
 
-  res = domain.integral('basis_ni,j stress_ij d:x' @ ns, degree=degree*2)
+  res = domain.integral('d(basis_ni, x_j) stress_ij d:x' @ ns, degree=degree*2)
   lhs = solver.solve_linear('lhs', res, constrain=cons)
 
   bezier = domain.sample('bezier', 5)

examples/finitestrain.py

@@ -46,7 +46,7 @@ def main(nelems:int, etype:str, btype:str, degree:int, poisson:float, angle:floa
   ns.ubasis = domain.basis(btype, degree=degree).vector(2)
   ns.u_i = 'ubasis_ki ?lhs_k'
   ns.X_i = 'x_i + u_i'
-  ns.strain_ij = '.5 (u_i,j + u_j,i)'
+  ns.strain_ij = '.5 (d(u_i, x_j) + d(u_j, x_i))'
   ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
 
   sqr = domain.boundary['left'].integral('u_k u_k d:x' @ ns, degree=degree*2)
@@ -58,7 +58,7 @@ def main(nelems:int, etype:str, btype:str, degree:int, poisson:float, angle:floa
   X, energy = bezier.eval(['X_i', 'energy'] @ ns, lhs=lhs0)
   export.triplot('linear.png', X, energy, tri=bezier.tri, hull=bezier.hull)
 
-  ns.strain_ij = '.5 (u_i,j + u_j,i + u_k,i u_k,j)'
+  ns.strain_ij = '.5 (d(u_i, x_j) + d(u_j, x_i) + d(u_k, x_i) d(u_k, x_j))'
   ns.energy = 'lmbda strain_ii strain_jj + 2 mu strain_ij strain_ij'
 
   energy = domain.integral('energy d:x' @ ns, degree=degree*2)

examples/laplace.py

@@ -61,13 +61,13 @@ def main(nelems:int, etype:str, btype:str, degree:int):
   # We are now ready to implement the Laplace equation. In weak form, the
   # solution is a scalar field :math:`u` for which:
   #
-  # .. math:: ∀ v: ∫_Ω v_{,k} u_{,k} - ∫_{Γ_n} v f = 0.
+  # .. math:: ∀ v: ∫_Ω \frac{dv}{dx_i} \frac{du}{dx_i} - ∫_{Γ_n} v f = 0.
   #
   # By linearity the test function :math:`v` can be replaced by the basis that
   # spans its space. The result is an integral ``res`` that evaluates to a
   # vector matching the size of the function space.
 
-  res = domain.integral('basis_n,i u_,i d:x' @ ns, degree=degree*2)
+  res = domain.integral('d(basis_n, x_i) d(u, x_i) d:x' @ ns, degree=degree*2)
   res -= domain.boundary['right'].integral('basis_n cos(1) cosh(x_1) d:x' @ ns, degree=degree*2)
 
   # The Dirichlet constraints are set by finding the coefficients that minimize

examples/platewithhole-nurbs.py

@@ -59,7 +59,7 @@ def main(nrefine:int, traction:float, radius:float, poisson:float):
   ns.ubasis = nurbsbasis.vector(2)
   ns.u_i = 'ubasis_ni ?lhs_n'
   ns.X_i = 'x_i + u_i'
-  ns.strain_ij = '(u_i,j + u_j,i) / 2'
+  ns.strain_ij = '(d(u_i, x_j) + d(u_j, x_i)) / 2'
   ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
   ns.r2 = 'x_k x_k'
   ns.R2 = radius**2 / ns.r2
@@ -74,7 +74,7 @@ def main(nrefine:int, traction:float, radius:float, poisson:float):
   cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
 
   # construct residual
-  res = domain.integral('ubasis_ni,j stress_ij d:x' @ ns, degree=9)
+  res = domain.integral('d(ubasis_ni, x_j) stress_ij d:x' @ ns, degree=9)
 
   # solve system
   lhs = solver.solve_linear('lhs', res, constrain=cons)
@@ -85,7 +85,7 @@ def main(nrefine:int, traction:float, radius:float, poisson:float):
   export.triplot('stressxx.png', X, stressxx, tri=bezier.tri, hull=bezier.hull, clim=(numpy.nanmin(stressxx), numpy.nanmax(stressxx)))
 
   # evaluate error
-  err = domain.integral('<du_k du_k, du_i,j du_i,j>_n d:x' @ ns, degree=9).eval(lhs=lhs)**.5
+  err = domain.integral('<du_k du_k, d(du_i, x_j) d(du_i, x_j)>_n d:x' @ ns, degree=9).eval(lhs=lhs)**.5
   treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
 
   return err, cons, lhs

examples/platewithhole.py

@@ -45,7 +45,7 @@ def main(nelems:int, etype:str, btype:str, degree:int, traction:float, maxrefine
   ns.ubasis = domain.basis(btype, degree=degree).vector(2)
   ns.u_i = 'ubasis_ni ?lhs_n'
   ns.X_i = 'x_i + u_i'
-  ns.strain_ij = '(u_i,j + u_j,i) / 2'
+  ns.strain_ij = '(d(u_i, x_j) + d(u_j, x_i)) / 2'
   ns.stress_ij = 'lmbda strain_kk δ_ij + 2 mu strain_ij'
   ns.r2 = 'x_k x_k'
   ns.R2 = radius**2 / ns.r2
@@ -59,14 +59,14 @@ def main(nelems:int, etype:str, btype:str, degree:int, traction:float, maxrefine
   sqr = domain.boundary['top,right'].integral('du_k du_k d:x' @ ns, degree=20)
   cons = solver.optimize('lhs', sqr, droptol=1e-15, constrain=cons)
 
-  res = domain.integral('ubasis_ni,j stress_ij d:x' @ ns, degree=degree*2)
+  res = domain.integral('d(ubasis_ni, x_j) stress_ij d:x' @ ns, degree=degree*2)
   lhs = solver.solve_linear('lhs', res, constrain=cons)
 
   bezier = domain.sample('bezier', 5)
   X, stressxx = bezier.eval(['X_i', 'stress_00'] @ ns, lhs=lhs)
   export.triplot('stressxx.png', X, stressxx, tri=bezier.tri, hull=bezier.hull)
 
-  err = domain.integral('<du_k du_k, du_i,j du_i,j>_n d:x' @ ns, degree=max(degree,3)*2).eval(lhs=lhs)**.5
+  err = domain.integral('<du_k du_k, d(du_i, x_j) d(du_i, x_j)>_n d:x' @ ns, degree=max(degree,3)*2).eval(lhs=lhs)**.5
   treelog.user('errors: L2={:.2e}, H1={:.2e}'.format(*err))
 
   return err, cons, lhs