Hi all. Consider the trivial code below that solves the tutorial example of the Poisson equation

```
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import FunctionSpace
from dolfinx import fem
import numpy as np
import ufl
import dolfinx
from petsc4py.PETSc import ScalarType
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
points = domain.geometry.x
V = FunctionSpace(domain, ("CG", 1))
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
bc = fem.dirichletbc(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType(-6))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
n=len(points)
xx=np.zeros(n)
yy=np.zeros(n)
zz=np.zeros(n)
for i in range(n):
xx[i]=points[i][0]
yy[i]=points[i][1]
zz[i]=points[i][2]
u=uh.x.array
```

After the calculation we basically have a matrix x(i),y(i),u(i), i=0,n-1 data points, where x(i) nad y(i) are the mesh points and u(i) is the solution on each point.

Consider now the gradient of u (or the partial derivative) on each mesh points. What is the easiest way to do this in ufl? I guess du=ufl.grad(u) but how can I save the result in the mesh points? Familiar method from the OG FEniCS do not work.

Thanks in advance