Here is an example of how to get the midpoints, and the values of u at the midpoints:
from mpi4py import MPI
import dolfinx
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.functionspace(mesh, ("CG", 2))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0] + x[1] ** 2)
Q = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.geometry.dim,)))
q = dolfinx.fem.Function(Q)
q.interpolate(lambda x: x[: mesh.geometry.dim])
num_cells_local = (
mesh.topology.index_map(mesh.topology.dim).size_local
+ mesh.topology.index_map(mesh.topology.dim).num_ghosts
)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
exterior_facet_indices = dolfinx.mesh.exterior_facet_indices(mesh.topology)
num_facets_local = mesh.topology.index_map(mesh.topology.dim - 1).size_local
interior_facet_indices = np.arange(num_facets_local)
interior_facet_indices = np.delete(interior_facet_indices, exterior_facet_indices)
import basix
ref_cell_geometry = basix.cell.geometry(mesh.basix_cell())
midpoint = np.sum(ref_cell_geometry, axis=0) / ref_cell_geometry.shape[0]
u_at_midpoint = dolfinx.fem.Expression(u, midpoint.reshape(-1, len(midpoint)))
midpoint_values = u_at_midpoint.eval(mesh, np.arange(num_cells_local, dtype=np.int32))
f_to_c = mesh.topology.connectivity(mesh.topology.dim - 1, mesh.topology.dim)
q_r = q.x.array.reshape(-1, mesh.geometry.dim)
for interior_facet in interior_facet_indices:
cells = f_to_c.links(interior_facet)
u_values = midpoint_values[cells]
q_values = q_r[cells]
assert np.isclose(u_values[0][0], q_values[0][0] + q_values[0][1] ** 2)
assert np.isclose(u_values[1][0], q_values[1][0] + q_values[1][1] ** 2)
This can now easily be extended to loop through each cell, then finding the facets of that cell, finding the neighbour, getting the midpoints to the adjacent cells and evaluate u in one cell minus u in another.