Measure function value at point

General
1

I have a fairly simple question: I want to measure the value of my function (u) at a certain point, let’s say (0.5, 0.5). How do I approach this? There are two cases, either the mesh intersects the probe point, or it does not, and I would like the solution to be compatible with both cases.
I have tried to use compute_closest_entity, but this gives me indices that are not in u.

Minimum working example (based on Poisson tutorial at https://github.com/jorgensd/dolfinx-tutorial/blob/main/chapter1/fundamentals_code.py):

from mpi4py import MPI
from dolfinx import mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

from dolfinx.fem import FunctionSpace
V = FunctionSpace(domain, ("CG", 1))

from dolfinx import fem
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)

import numpy
# Create facet to cell connectivity required to determine boundary facets
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)

import ufl
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

from petsc4py.PETSc import ScalarType
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()

import pyvista
print(pyvista.global_theme.jupyter_backend)

# + vscode={"languageId": "python"}
from dolfinx import plot
pyvista.start_xvfb()
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(V)

u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh.x.array.real
u_grid.set_active_scalars("u")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()

# ```

What should I add to this MWE to measure the value of uh at point (0.5, 0.5)?

1 Like
2

Could you please provide the attempted solution, as you have not added any code about evaluating at points to this example?

Please also note that it would be sufficient to use the following lines to illustrate your issue:

from mpi4py import MPI
from dolfinx import mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

from dolfinx.fem import FunctionSpace
V = FunctionSpace(domain, ("CG", 1))

from dolfinx import fem
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)

# Add code here to evaluate uD at a given point (say 0.5, 0.5)
# ...

Have you also considered: Implementation — FEniCSx tutorial ?

3

Thank you for the quick response! Would you say using compute_closest_entity is the recommended method for doing this?

This is a MWE with my attempt at using compute_closest_entity.

from mpi4py import MPI
from dolfinx import mesh
import numpy as np
from dolfinx.geometry import (BoundingBoxTree, compute_closest_entity,
                              create_midpoint_tree)
from pdb import set_trace


domain = mesh.create_unit_square(MPI.COMM_WORLD, 80, 80, mesh.CellType.quadrilateral)

from dolfinx.fem import FunctionSpace
V = FunctionSpace(domain, ("CG", 1))

from dolfinx import fem
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0])
dim = domain.topology.dim

# from https://github.com/FEniCS/dolfinx/blob/a5b84c10e8ccc6d236f780d3a0486eb634d048c0/python/test/unit/geometry/test_bounding_box_tree.py#L277
points = np.array([1.0, -1.0])
domain.topology.create_entities(dim)
# set_trace()
tree = BoundingBoxTree(domain, dim)
num_entities_local = domain.topology.index_map(dim).size_local + domain.topology.index_map(dim).num_ghosts
entities = np.arange(num_entities_local, dtype=np.int32)
midpoint_tree = create_midpoint_tree(domain, dim, entities)

closest_entities = compute_closest_entity(tree, midpoint_tree, domain, points)

# get value at closest_entity index
value = uD.x.array[closest_entities]
print(f" value at {points}: {value}")

The expected value at the evaluated point is x[0]+1.0=2.0, but the output is:
value at [ 1. -1.]: [1.0125]

I can imagine there are simpler methods to measure values at a point. The tutorial shows many ways to plot the values in pyvista, which works for me, but I would like to create a graph of the change of u over time for some points in my simulation.
Thanks again for the help.

4

You are mixing up entities and degrees of freedom.

Entities of dimension d are the components of a mesh of that dimension.

  • d=0 is points
  • d=1 is intervals
  • d=2 is triangles/quadrilaterals
  • d=3 is tetrahedra/hexahedra

The mesh.topology.dim tells you what dimension the cells of the mesh have.
So what you are currently doing is:

  1. Finding the cell c_i that is closest to ([1, -1]) (which is not in your mesh btw).
  2. Evaluating the c_ith dof of your function (which is in a first order Lagrange space, which has no connection to the cell index).

The tutorial also shows how to plot u over a line. i.e. you can reduce that code to plotting over a single point in a time loop

image
image1499×1839 265 KB

image
image1499×1269 266 KB

1 Like
5

Thank you, simplifying the line method to a point seems to be exactly what I need.