Finding facet midpoint on mesh using bounding box tree

mesh
1

Hi,

consider this MWE using Fenics 2017.2.0:

from dolfin import *
import numpy as np
from mshr import *

__UINT32_MAX__ = np.iinfo('uint32').max

domain = Sphere(Point(0,0,0), 1)
mesh = generate_mesh(domain, 32)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(40)

tree = mesh.bounding_box_tree()

facets = SubsetIterator(boundaries, 40)
count = 0
for facet in facets:
    centre = facet.midpoint()
    collision = tree.compute_first_entity_collision(centre)
    if collision != -1 and collision != __UINT32_MAX__:
        count += 1

print(count)

I get the following output

Generating mesh with CGAL 3D mesh generator
Iterating over subset, found 153991 entities out of 153991.
151412

which basically tells me that the number stored in the variable count is less than the number of facet entities in the mesh. I need to find all facet midpoints in the mesh, why does this method not seem to work? How can I fix it?

2

I try to understand your problem but it is not clear to me why you still need to do the collision at the first place if you can already access the facet midpoints via facet.midpoint(). On the other hand, there might be some numerical issues with the compute_first_entity_collision so that some centre points are not found.