I’m having trouble with the dolfin.cpp.mesh.Cell.collides() function, as it seems it has a huge likelyhood of false negatives, but it seems rather “random”.
I wanted to use it to be able to locate a point inside the mesh, and later interpolate the value of fields at the point using the values of the cell containing the point (and I cannot use things like project() because it doesn’t work properly in a multithread environment)
Below is an example showing the issue. It is supposed to be ran with fenics’s latest stable version, mpi and python 3:
mpirun -n X python3 test.py
You should see some lines saying that some points were not seen inside the mesh, and if you change the number of threads it runs on, or even if you rerun it, you will probably have a different output.
Is there a way to achieve what I want without re-implementing the “point inside a segment/triangle/tetrahedron” methods myself ?
from dolfin import *
comm = MPI.comm_world
rank = MPI.rank(comm)
size_mpi= MPI.size(comm)
p1 = Point(-100,-100,0)
p2 = Point(100,100,0)
m = RectangleMesh(p1, p2, 200,200)
v = FunctionSpace(m, "P", 1)
d = v.dofmap()
my_first, my_last = d.ownership_range()
#Creation of a list of points that are geometrically inside the mesh
pts = [Point(x,x,0) for x in range(100)]
pt_collision = [0 for x in pts]
pt_owner = [-1 for x in pts]
print("rank " + str(rank) + " " + str(my_first) + " " + str(my_last))
for cell in cells(m):
for i, p in enumerate(pts):
# Test if the point in question in inside the cell
collision = cell.collides(p)
collision_number = 0
if collision:
# Test if the cell is completly owned by the current thread
is_fully_owned = True
for vert in vertices(cell):
ind = vert.index()
global_index = d.local_to_global_index(ind)
if not (my_first <= global_index < my_last):
is_fully_owned = False
if is_fully_owned and pt_collision[i] == 0:
pt_collision[i] =1
pt_owner[i] = rank
print("in rank " + str(rank))
print("pt of coor " + str(p[0])+ " " +str(p[1]))
print("inside : ")
for vert in vertices(cell):
print(str(vert.point()[0]) + " " +str(vert.point()[1]) + " with global index " + str(d.local_to_global_index(vert.index())) + " local index " + str(vert.index()))
print()
print()
for point_index, point_collisions in enumerate(pt_collision):
owner_index = MPI.max(comm, pt_owner[point_index])
if owner_index == -1 and rank == 0:
print("a point was not found")
for r in range(size_mpi):
if rank == r and r!= owner_index and pt_owner[point_index] == r:
pt_collision[point_index] = 0
print("a point was seen more than once")
Hi, have you tried adapting the implementation of Probe class of fenicstools library? A probe at given point works only with dofs local to the cell where the point is located; this sounds like what your are after. In addition the implementation for many probes is \mathcal{O}(\#points) wheras yours is \mathcal{O}(\#cells).
I just tried. fenicstools’ Probe use mesh.bounding_box_tree().compute_first_entity_collision() to get the cell, so I tried that, and still i get some points that are not found in the mesh.
Besides, I would greatly expect compute_first_entity_collision() to internally use cell.collides(), though I did not check the internal code to verify this.
I did not install fenicstools, because I couldn’t get it to work (imports would fail, blaming missing dependencies).
Here is the modified code I came up with.
from dolfin import *
comm = MPI.comm_world
rank = MPI.rank(comm)
size_mpi= MPI.size(comm)
p1 = Point(-100,-100,0)
p2 = Point(100,100,0)
m = RectangleMesh(p1, p2, 200,200)
v = FunctionSpace(m, "P", 1)
d = v.dofmap()
my_first, my_last = d.ownership_range()
pts = [Point(x,x,0) for x in range(100)]
pt_collision = [0 for x in pts]
pt_owner = [-1 for x in pts]
print("rank " + str(rank) + " " + str(my_first) + " " + str(my_last))
#for cell in cells(m):
for i, p in enumerate(pts):
#collision = cell.collides(p)
#collision_number = 0
cell_index = m.bounding_box_tree().compute_first_entity_collision(p)
if cell_index < 2**30:
is_fully_owned = True
for vert in m.cells()[cell_index]:
#ind = vert.index()
ind = vert
global_index = d.local_to_global_index(ind)
if not (my_first <= global_index < my_last):
is_fully_owned = False
if is_fully_owned and pt_collision[i] == 0:
pt_collision[i] =1
pt_owner[i] = rank
#print("in rank " + str(rank))
#print("pt of coor " + str(p[0])+ " " +str(p[1]))
#print("inside : ")
#for vert in vertices(cell):
#for vert in m.cells()[cell_index]:
# print(str(vert.point()[0]) + " " +str(vert.point()[1]) + " with global index " + str(d.local_to_global_index(vert.index())) + " local index " + str(vert.index()))
#print()
#print()
for point_index, point_collisions in enumerate(pt_collision):
owner_index = MPI.max(comm, pt_owner[point_index])
if owner_index == -1 and rank == 0:
print("a point was not found")
for r in range(size_mpi):
if rank == r and r!= owner_index and pt_owner[point_index] == r:
pt_collision[point_index] = 0
print("a point was seen more than once")
PS: the complexity of the algorithm is not a critical issue here, i’ll be doing this only once during hour-long simulations, so i’m ok with not-so optimized solutions
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
I have tried your code and want to apply it to vector function. It produced wrong results.
I just change the function space, expression of f, size of A , and last several lines. I can not find where is the mistake.
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 = VectorFunctionSpace(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 = VectorFunctionSpace(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()*2)
# 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)
value1 = coefs.dot(A[0:6])
value2 = coefs.dot(A[6:12])
x, y, _ = pt.array()
print("error_1: ", value1-2*x)
print("error_2: ", value2-3*y)
I have solved the problem. Because “evaluate_basis_all” return a column first matrix, transpose it to get the right answer. Here is my code, just substitute several lines.