The points have the wrong shape (single point, the function expects multiple points):
from mpi4py import MPI
import numpy as np
from dolfinx import mesh
from dolfinx import fem
from dolfinx.fem import FunctionSpace
from dolfinx import default_scalar_type
import ufl
from dolfinx.fem.petsc import LinearProblem
from dolfinx import geometry
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("Lagrange", 1))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type(-6))
# creation of the bilinear form
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(
a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
uh = problem.solve()
points = np.array([[0.5, 0.5, 0.0]])
bb_tree = geometry.bb_tree(domain, domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions_points(bb_tree, points)
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points)
for i, point in enumerate(points):
if len(colliding_cells.links(i)) > 0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
u_val = uh.eval(points_on_proc, cells)
print(u_val)
However, I would rather go for something along the lines of:
i.e.
from mpi4py import MPI
import numpy as np
from dolfinx import mesh
from dolfinx import fem, cpp, default_scalar_type
from dolfinx.fem import FunctionSpace
from dolfinx import default_scalar_type
import ufl
from dolfinx.fem.petsc import LinearProblem
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("Lagrange", 2))
def f(x):
return x[0] ** 2 + 3 * x[1] * x[1]
def gradf(x):
return np.asarray((2 * x[0], 6 * x[1]))
uh = fem.Function(V)
uh.interpolate(f)
if MPI.COMM_WORLD.rank == 0:
points = np.array(
[[0.5, 0.5, 0.0], [0.25, 0.76, 0], [0.64, 0.33, 0]],
dtype=domain.geometry.x.dtype,
)
else:
points = np.empty((0, 3), dtype=domain.geometry.x.dtype)
def evaluate_expression(mesh, expr, points):
# Determine what process owns a point and what cells it lies within
_, _, owning_points, cells = cpp.geometry.determine_point_ownership(
mesh._cpp_object, points, 1e-6
)
owning_points = np.asarray(owning_points).reshape(-1, 3)
# Pull owning points back to reference cell
mesh_nodes = mesh.geometry.x
cmap = mesh.geometry.cmaps[0]
ref_x = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh.geometry.x.dtype)
for i, (point, cell) in enumerate(zip(owning_points, cells)):
geom_dofs = mesh.geometry.dofmap[cell]
ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])
# Create expression evaluating a trial function (i.e. just the basis function)
data_size = np.prod(expr.ufl_shape)
if len(cells) > 0:
# NOTE: Expression lives on only this communicator rank
expr = fem.Expression(expr, ref_x, comm=MPI.COMM_SELF)
values = expr.eval(mesh, np.asarray(cells, dtype=np.int32))
# strip basis_values per cell
basis_values = np.empty((len(cells), data_size), dtype=default_scalar_type)
for i in range(len(cells)):
basis_values[i] = values[i, i * data_size : (i + 1) * data_size]
else:
basis_values = np.zeros((0, data_size), dtype=default_scalar_type)
return cells, basis_values
print(
MPI.COMM_WORLD.rank,
evaluate_expression(domain, ufl.grad(uh), points)[1],
flush=True,
)
MPI.COMM_WORLD.Barrier()
if MPI.COMM_WORLD.rank == 0:
print(gradf(points.T).T)
producing
0 [[1. 3. ]
[0.5 4.56]]
1 [[1.28 1.98]]
[[1. 3. ]
[0.5 4.56]
[1.28 1.98]]