dolfin.cpp.mesh.Cell.collides()

mesh
4

Okay but the process does not need to own all of the cell’s dofs for interpolation to work. This is shown below using ideas of the fenicstools.Probe

from dolfin import *
import numpy as np

comm = MPI.comm_world
rank = MPI.rank(comm)
size_mpi = MPI.size(comm)

mesh = UnitSquareMesh(128, 128)
v = FunctionSpace(mesh, 'CG', 1)

pts = np.r_[np.random.rand(200, 2), np.array([[2, 2]])] if rank == 0 else None
# Everybody get same (root) pts
pts = comm.bcast(pts)
    
pts = list(map(Point, pts))
cells_of_points = [-1]*len(pts)

bbox_tree, ncells = mesh.bounding_box_tree(), mesh.num_cells()
for i, p in enumerate(pts):
    cell_index = bbox_tree.compute_first_collision(p)
    if cell_index < ncells:
        cells_of_points[i] = cell_index
# Who found what
gcells_of_points = comm.allreduce([cells_of_points])

is_lost = lambda collection: size_mpi == sum(1 for i in collection if i == -1)

lost = list(map(is_lost, zip(*gcells_of_points)))
assert sum(lost) == 1

is_local = lambda collection: collection[rank] != -1
local = list(map(is_local, zip(*gcells_of_points)))

# We can interpolate even if cell does not own all dofs on it
V = FunctionSpace(mesh, 'CG', 1)
f = interpolate(Expression(('2*x[0]-3*x[1]'), degree=1), V)
vec = as_backend_type(f.vector()).vec()

element = V.dolfin_element()
dm = V.dofmap()

A = np.zeros(element.space_dimension())
# Get the value of f(pt) by restricting coef vector to cell and dotting
# the local vector with matrix represening cell basis function values at pt
# This is what f.interpolate does but without the search for cells
for i, (pt, is_local) in enumerate(zip(pts, local)):
    if is_local:
        cell = cells_of_points[i]
        # Global  = local + offset
        cell_dofs = dm.cell_dofs(cell) + dm.ownership_range()[0]
        # Basis foo values aat pt
        cell = Cell(mesh, cell)
        vertex_x, orientation = cell.get_vertex_coordinates(), cell.orientation()
        A[:] = element.evaluate_basis_all(pt.array(), vertex_x, orientation)
        # Element sum 
        coefs = vec.getValues(cell_dofs)
        value = coefs.dot(A)

        x, y, _ = pt.array()
        # Our reference is sum of pt as f = 2x-3y
        assert abs(value - (2*x-3*y)) < 1E-14
1 Like