Suppose I’m solving a BVP
L(u) = f, \quad \mathrm{in} \; \Omega \subset \mathbb{R}^2
with
u = g_D, \quad \mathrm{on} \; \partial\Omega_D
and I want to evaluate u, \nabla u at some \mathbf{x} \in \Omega. Is there an easy way to do this in FEniCS or do I have to manually get the values at DoFs and perform the interpolation between them?
This depends on if you are using fenics or FEniCSx.
See for instance: Evaluation of UFL expression at mesh points - #6 by dokken and subsequent points.
Great, this handles the evaluation of u. Now, is there any easy way to evaluate \nabla u at an arbitrary point \mathbf{x} \in \Omega?
I attempted to copy the procedure shown here https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain to evaluate a finite element function by appending the following lines
point = np.array([0.5,0.5])
bb_tree = geometry.bb_tree(domain, domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions_points(bb_tree, point.T)
colliding_cells = geometry.compute_colliding_cells(domain,cell_candidates,point.T)
for i, point in enumerate(point.T):
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)
To this tutorial https://jsdokken.com/dolfinx-tutorial/chapter1/fundamentals_code.html.
However, when I do this, the following error message is printed to console,
AttributeError: 'numpy.ndarray' object has no attribute 'links'
Referring to
if len(colliding_cells.links(i)) > 0:
In my approach, `colliding_cells’ is saved as a numpy array. Is this because the mesh was not generated with gmsh? If so, must one use gmsh in order to evaluate FE functions?
No requirement for using gmsh for this approach to work.
Please state what version of DOLFINx that you are using.
The tutorial is based on v0.7.3
https://jsdokken.com/dolfinx-tutorial/index.html
Any older versions can be found at:
I’m using version 0.7.3.
Please provide a minimal reproducible example I could run to debug your issue.
Here you are
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
import pyvista
from dolfinx import geometry
from dolfinx import plot
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()
point = np.array([0.5,0.5,0.])
bb_tree = geometry.bb_tree(domain, domain.topology.dim)
cells = []
points_on_proc = []
cell_candidates = geometry.compute_collisions_points(bb_tree, point.T)
colliding_cells = geometry.compute_colliding_cells(domain,cell_candidates,point.T)
for i, point in enumerate(point.T):
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)
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]]
Excellent, thank you! Much appreciated.