[MOAB-dev] Entity picking

[Jason Kraftcheck]

Roman Putanowicz wrote:
> Hi Everybody,
> 
> I would like to ask you for advice.
> 
> The task is: given a point position find mesh cell containing
> that point (or cells in case of boundary point). I understand
> that there is a delicacy in the above due to the inexact nature 
> of geometric queries in double precision arithmetic, but approximate
> solution based on some arbitrary tolerance will do.
> 
> The question is how the above could be solved using MOAB?
> 
> Looking at the class list I noticed MBOrientedBoxTreeTool,
> MBAdaptiveKDTree, MBBSPTree classes which probably could be used
> for that purpose. (MBCoupler??). The answer depends for sure on 
> the context -- at the moment I need to pick elements in medium 
> to small size, unstructured polygonal meshes in 2D.
> 

MBAdaptiveKDTree is likely the best choice for this.  However, it currently
doesn't work for general polygons (it will work for triangles and
quadrilaterals).  The box_elem_overlap function in MBGeomUtil.cpp needs to
handle the general polygon case.


> I will be grateful for any advice or code snippet if somebody has any.
> 

// build tree

MBRange polygons;
MBErrorCode rval;
rval = moab_instance->get_entities_by_type( 0, MBPOLYGON, polygons );
if (MB_SUCCESS != rval) ...

MBAdaptiveKDTreeTool tool( moab_instance );
MBEntityHandle tree_root;

rval = tool.build_tree( polygons, tree_root );
if (MB_SUCCESS != rval) ...

// get some polygons near a point P

MBEntityHandle poly_set;
rval = tool.leaf_containing_point( tree_root, P, poly_set );
//
// or perhaps:
//
// std::vector<MBEntityHandle> poly_sets;
// rval = tool.leaves_within_distance( tree_root, P, epsilon, poly_sets );
//
if (MB_SUCESS != rval) ...

polygons.clear();
rval = moab_instance->get_entities_by_handle( poly_set, polygons );

// test each polygon to see if it contains the point

std::vector<double> coords;
for (MBRange::iterator i = polygons.begin(); i != polygons.end(); ++i)
{
  const MBEntityHandle* conn;
  int len;
  rval = moab_instance->get_connectivity( *i, conn, len );
  if (MB_SUCCESS != rval) ...

  coords.resize( 3*len );
  rval = moab_instance->coords( conn, len, &coords[0] );

  ...
}