Function and expression evaluation#

In this section, we will look at how we can extract data once we have solved our variational problem. To make this section concise, we will use functions with expressions that are already defined (not through a PDE), but through interpolation.

We start by creating a 3D mesh

from mpi4py import MPI
import dolfinx
import numpy as np

N = 2
mesh = dolfinx.mesh.create_box(
    MPI.COMM_WORLD,
    [np.array([0, 0, 0]), np.array([2, 1.3, 0.8])],
    [N, N, N],
    dolfinx.mesh.CellType.tetrahedron,
    ghost_mode=dolfinx.mesh.GhostMode.shared_facet,
)

We start by considering a scalar function in a discontinuous space

V = dolfinx.fem.functionspace(mesh, ("DG", 2))
u = dolfinx.fem.Function(V)

Interpolation on a subset of cells#

We would like to interpolate the function

\[\begin{split} \begin{align} u = \begin{cases} x[0]^2 & \text{if } x[0] < 1\\ x[1] & \text{otherwise} \end{cases} \end{align} \end{split}\]

We start by locating the cells that satisfies this condition. In Locating a subset of entities on a boundary we learnt how to use dolfinx.mesh.locate_entities_boundary to locate entites. In this section we will use dolfinx.mesh.locate_entities which does the same, but for all entities in the mesh.

tdim = mesh.topology.dim
left_cells = dolfinx.mesh.locate_entities(mesh, tdim, lambda x: x[0] <= 1 + 1e-14)
right_cells = dolfinx.mesh.locate_entities(mesh, tdim, lambda x: x[0] >= 1 - 1e-14)

We can now interpolate the function onto each of these subsets

u.interpolate(lambda x: x[0] ** 2, cells0=left_cells)
u.interpolate(lambda x: x[1], cells0=right_cells)

Whenever we interpolate on sub-sets of cells, we need to scatter forward the values

u.x.scatter_forward()

Evalation at a point#

We want to evaluate the function at a point in the mesh that does not align with the nodes of the mesh. We do this through a sequence of steps:

  1. As the mesh will be distributed on multiple processes, we have to determine which processes that has the point.

  2. Then, for the processes that has the point, we need to figure out what cell the point is in.

  3. Finally, we push this point back to the reference element, so that we can evaluate the basis functions, combine them with the coefficients on the given cell and push them forward to the physical space.

Step 1: Determine which processes that has the cell#

As looping through all the cells, and compute exact collisions is expensive, we accelerate the search by using an axis-aligned bounding box tree. This is a tree that recursively divides the mesh into smaller and smaller boxes, such that we can quickly search through the tree to find the cells that might contain the point of interest.

bb_tree = dolfinx.geometry.bb_tree(mesh, mesh.topology.dim)

Note

Bounding boxes of other entities As seen below, we send in the topological dimension of the entities that we are interested in. This means that you can make a bounding box tree for facets, edges or vertices as well if it is needed.