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