Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
SpatialLocator.cpp
Go to the documentation of this file.
2 #include "moab/Interface.hpp"
3 #include "moab/ElemEvaluator.hpp"
5 #include "moab/BVHTree.hpp"
6 
7 // include ScdInterface for box partitioning
8 #include "moab/ScdInterface.hpp"
9 
10 #ifdef MOAB_HAVE_MPI
11 #include "moab/ParallelComm.hpp"
12 #endif
13 
14 namespace moab
15 {
16 static bool debug = false;
17 
19  : mbImpl( impl ), myElems( elems ), myDim( -1 ), myTree( tree ), elemEval( eval ), iCreatedTree( false ),
20  timerInitialized( false )
21 {
22  create_tree();
23 
24  if( !elems.empty() )
25  {
28  if( MB_SUCCESS != rval ) throw rval;
29 
30  rval = myTree->get_bounding_box( localBox );
31  if( MB_SUCCESS != rval ) throw rval;
32  }
33 
34  regNums[0] = regNums[1] = regNums[2] = 0;
35 }
36 
38 {
39  if( myTree ) return;
40 
42  // create a kdtree if only vertices
43  myTree = new AdaptiveKDTree( mbImpl );
44  else
45  // otherwise a BVHtree, since it performs better for elements
46  myTree = new BVHTree( mbImpl );
47 
48  iCreatedTree = true;
49 }
50 
52 {
53  if( elems.empty() ||
55  return MB_FAILURE;
56 
57  myDim = mbImpl->dimension_from_handle( *elems.begin() );
58  myElems = elems;
59 
61  return rval;
62 }
63 
64 #ifdef MOAB_HAVE_MPI
65 ErrorCode SpatialLocator::initialize_intermediate_partition( ParallelComm* pc )
66 {
67  if( !pc ) return MB_FAILURE;
68 
69  BoundBox gbox;
70 
71  // step 2
72  // get the global bounding box
73  double sendbuffer[6];
74  double rcvbuffer[6];
75 
76  localBox.get( sendbuffer ); // fill sendbuffer with local box, max values in [0:2] min values in [3:5]
77  sendbuffer[0] *= -1;
78  sendbuffer[1] *= -1; // negate Xmin,Ymin,Zmin to get their minimum using MPI_MAX
79  sendbuffer[2] *= -1; // to avoid calling MPI_Allreduce again with MPI_MIN
80 
81  int mpi_err = MPI_Allreduce( sendbuffer, rcvbuffer, 6, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
82  if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
83 
84  rcvbuffer[0] *= -1;
85  rcvbuffer[1] *= -1; // negate Xmin,Ymin,Zmin again to get original values
86  rcvbuffer[2] *= -1;
87 
88  globalBox.update_max( &rcvbuffer[3] ); // saving values in globalBox
89  globalBox.update_min( &rcvbuffer[0] );
90 
91  // compute the alternate decomposition; use ScdInterface::compute_partition_sqijk for this
92  ScdParData spd;
94  spd.gPeriodic[0] = spd.gPeriodic[1] = spd.gPeriodic[2] = 0;
95  double lg = log10( ( localBox.bMax - localBox.bMin ).length() );
96  double mfactor = pow( 10.0, 6 - lg );
97  int ldims[6], lper[3];
98  double dgijk[6];
99  localBox.get( dgijk );
100  for( int i = 0; i < 6; i++ )
101  spd.gDims[i] = dgijk[i] * mfactor;
102  ErrorCode rval = ScdInterface::compute_partition( pc->size(), pc->rank(), spd, ldims, lper, regNums );
103  if( MB_SUCCESS != rval ) return rval;
104  // all we're really interested in is regNums[i], #procs in each direction
105 
106  for( int i = 0; i < 3; i++ )
107  regDeltaXYZ[i] = ( globalBox.bMax[i] - globalBox.bMin[i] ) / double( regNums[i] ); // size of each region
108 
109  return MB_SUCCESS;
110 }
111 
112 // this function sets up the TupleList TLreg_o containing the registration messages
113 // and sends it
114 ErrorCode SpatialLocator::register_src_with_intermediate_procs( ParallelComm* pc,
115  double abs_iter_tol,
116  TupleList& TLreg_o )
117 {
118  int corner_ijk[6];
119 
120  // step 3: compute ijks of local box corners in intermediate partition
121  // get corner ijk values for my box
122  ErrorCode rval = get_point_ijk( localBox.bMin - CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk );
123  if( MB_SUCCESS != rval ) return rval;
124  rval = get_point_ijk( localBox.bMax + CartVect( abs_iter_tol ), abs_iter_tol, corner_ijk + 3 );
125  if( MB_SUCCESS != rval ) return rval;
126 
127  // step 4
128  // set up TLreg_o
129  TLreg_o.initialize( 1, 0, 0, 6, 0 );
130  // TLreg_o (int destProc, real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax)
131 
132  int dest;
133  double boxtosend[6];
134 
135  localBox.get( boxtosend );
136 
137  // iterate over all regions overlapping with my bounding box using the computerd corner IDs
138  for( int k = corner_ijk[2]; k <= corner_ijk[5]; k++ )
139  {
140  for( int j = corner_ijk[1]; j <= corner_ijk[4]; j++ )
141  {
142  for( int i = corner_ijk[0]; i <= corner_ijk[3]; i++ )
143  {
144  dest = k * regNums[0] * regNums[1] + j * regNums[0] + i;
145  TLreg_o.push_back( &dest, NULL, NULL, boxtosend );
146  }
147  }
148  }
149 
150  // step 5
151  // send TLreg_o, receive TLrequests_i
152  if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLreg_o, 0 );
153 
154  // step 6
155  // Read registration requests from TLreg_o and add to list of procs to forward to
156  // get number of tuples sent to me
157 
158  // read tuples and fill processor list;
159  int NN = TLreg_o.get_n();
160  for( int i = 0; i < NN; i++ )
161  // TLreg_o is now TLrequests_i
162  srcProcBoxes[TLreg_o.vi_rd[i]] = BoundBox( TLreg_o.vr_rd + 6 * i );
163 
164  return MB_SUCCESS;
165 }
166 
167 ErrorCode SpatialLocator::par_locate_points( ParallelComm* /*pc*/,
168  Range& /*vertices*/,
169  const double /*rel_iter_tol*/,
170  const double /*abs_iter_tol*/,
171  const double /*inside_tol*/ )
172 {
174 }
175 
176 bool is_neg( int i )
177 {
178  return ( i == -1 );
179 }
180 
181 ErrorCode SpatialLocator::par_locate_points( ParallelComm* pc,
182  const double* pos,
183  int num_points,
184  const double rel_iter_tol,
185  const double abs_iter_tol,
186  const double inside_tol )
187 {
188  ErrorCode rval;
189  // TUpleList used for communication
190  TupleList TLreg_o; // TLregister_outbound
191  TupleList TLquery_o; // TLquery_outbound
192  TupleList TLforward_o; // TLforward_outbound
193  TupleList TLsearch_results_o; // TLsearch_results_outbound
194 
195  // initialize timer
197  timerInitialized = true;
198 
199  // steps 1-2 - initialize the alternative decomposition box from global box
200  rval = initialize_intermediate_partition( pc );
201  if( rval != MB_SUCCESS ) return rval;
202 
203  // steps 3-6 - set up TLreg_o, gs_transfer, gather registrations
204  rval = register_src_with_intermediate_procs( pc, abs_iter_tol, TLreg_o );
205  if( rval != MB_SUCCESS ) return rval;
206 
208 
209  // actual parallel point location using intermediate partition
210 
211  // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing
212  // pts to be located source_pts: TL(from_proc, tgt_index, src_index): results of source mesh
213  // proc point location, ready to send
214  // back to tgt procs; src_index of -1 indicates point not located (arguably not
215  // useful...)
216 
217  unsigned int my_rank = ( pc ? pc->proc_config().proc_rank() : 0 );
218 
219  // TLquery_o: Tuples sent to forwarder proc
220  // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
221 
222  // TLforw_req_i: Tuples to forward to corresponding procs (forwarding requests)
223  // TL (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
224 
225  TLquery_o.initialize( 3, 0, 0, 3, 0 );
226 
227  int iargs[3];
228 
229  for( int pnt = 0; pnt < 3 * num_points; pnt += 3 )
230  {
231  int forw_id =
232  proc_from_point( pos + pnt, abs_iter_tol ); // get ID of proc resonsible of the region the proc is in
233 
234  iargs[0] = forw_id; // toProc
235  iargs[1] = my_rank; // originalSourceProc
236  iargs[2] = pnt / 3; // targetIndex
237 
238  TLquery_o.push_back( iargs, NULL, NULL, const_cast< double* >( pos + pnt ) );
239  }
240 
241  // send point search queries to forwarders
242  if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLquery_o, 0 );
243 
245 
246  // now read forwarding requests and forward to corresponding procs
247  // TLquery_o is now TLforw_req_i
248 
249  // TLforward_o: query messages forwarded to corresponding procs
250  // TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
251 
252  TLforward_o.initialize( 3, 0, 0, 3, 0 );
253 
254  int NN = TLquery_o.get_n();
255 
256  for( int i = 0; i < NN; i++ )
257  {
258  iargs[1] = TLquery_o.vi_rd[3 * i + 1]; // get OriginalSourceProc
259  iargs[2] = TLquery_o.vi_rd[3 * i + 2]; // targetIndex
260  CartVect tmp_pnt( TLquery_o.vr_rd + 3 * i );
261 
262  // compare coordinates to list of bounding boxes
263  for( std::map< int, BoundBox >::iterator mit = srcProcBoxes.begin(); mit != srcProcBoxes.end(); ++mit )
264  {
265  if( ( *mit ).second.contains_point( tmp_pnt.array(), abs_iter_tol ) )
266  {
267  iargs[0] = ( *mit ).first;
268  TLforward_o.push_back( iargs, NULL, NULL, tmp_pnt.array() );
269  }
270  }
271  }
272 
274 
275  if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLforward_o, 0 );
276 
278 
279  // cache time here, because locate_points also calls elapsed functions and we want to account
280  // for tuple list initialization here
281  double tstart = myTimer.time_since_birth();
282 
283  // step 12
284  // now read Point Search requests
285  // TLforward_o is now TLsearch_req_i
286  // TLsearch_req_i: (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
287 
288  NN = TLforward_o.get_n();
289 
290  // TLsearch_results_o
291  // TL: (OriginalSourceProc, targetIndex, sourceIndex, U,V,W);
292  TLsearch_results_o.initialize( 3, 0, 0, 0, 0 );
293 
294  // step 13 is done in test_local_box
295 
296  std::vector< double > params( 3 * NN );
297  std::vector< int > is_inside( NN, 0 );
298  std::vector< EntityHandle > ents( NN, 0 );
299 
300  rval = locate_points( TLforward_o.vr_rd, TLforward_o.get_n(), &ents[0], &params[0], &is_inside[0], rel_iter_tol,
301  abs_iter_tol, inside_tol );
302  if( MB_SUCCESS != rval ) return rval;
303 
304  locTable.initialize( 1, 0, 1, 3, 0 );
306  for( int i = 0; i < NN; i++ )
307  {
308  if( is_inside[i] )
309  {
310  iargs[0] = TLforward_o.vi_rd[3 * i + 1];
311  iargs[1] = TLforward_o.vi_rd[3 * i + 2];
312  iargs[2] = locTable.get_n();
313  TLsearch_results_o.push_back( iargs, NULL, NULL, NULL );
314  Ulong ent_ulong = (Ulong)ents[i];
315  sint forward = (sint)TLforward_o.vi_rd[3 * i + 1];
316  locTable.push_back( &forward, NULL, &ent_ulong, &params[3 * i] );
317  }
318  }
320 
322  myTimer.time_elapsed(); // call this to reset last time called
323 
324  // step 14: send TLsearch_results_o and receive TLloc_i
325  if( pc ) pc->proc_config().crystal_router()->gs_transfer( 1, TLsearch_results_o, 0 );
326 
328 
329  // store proc/index tuples in parLocTable
330  parLocTable.initialize( 2, 0, 0, 0, num_points );
332  std::fill( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, -1 );
333 
334  for( unsigned int i = 0; i < TLsearch_results_o.get_n(); i++ )
335  {
336  int idx = TLsearch_results_o.vi_rd[3 * i + 1];
337  parLocTable.vi_wr[2 * idx] = TLsearch_results_o.vi_rd[3 * i];
338  parLocTable.vi_wr[2 * idx + 1] = TLsearch_results_o.vi_rd[3 * i + 2];
339  }
340 
341  if( debug )
342  {
343  int num_found =
344  num_points - 0.5 * std::count_if( parLocTable.vi_wr, parLocTable.vi_wr + 2 * num_points, is_neg );
345  std::cout << "Points found = " << num_found << "/" << num_points << " ("
346  << 100.0 * ( (double)num_found / num_points ) << "%)" << std::endl;
347  }
348 
350 
351  return MB_SUCCESS;
352 }
353 
354 #endif
355 
357  const double rel_iter_tol,
358  const double abs_iter_tol,
359  const double inside_tol )
360 {
361  bool i_initialized = false;
362  if( !timerInitialized )
363  {
365  timerInitialized = true;
366  i_initialized = true;
367  }
368 
369  assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
370  std::vector< double > pos( 3 * verts.size() );
371  ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
372  if( MB_SUCCESS != rval ) return rval;
373  rval = locate_points( &pos[0], verts.size(), rel_iter_tol, abs_iter_tol, inside_tol );
374  if( MB_SUCCESS != rval ) return rval;
375 
376  // only call this if I'm the top-level function, since it resets the last time called
378 
379  return MB_SUCCESS;
380 }
381 
383  int num_points,
384  const double rel_iter_tol,
385  const double abs_iter_tol,
386  const double inside_tol )
387 {
388  bool i_initialized = false;
389  if( !timerInitialized )
390  {
392  timerInitialized = true;
393  i_initialized = true;
394  }
395  // initialize to tuple structure (p_ui, hs_ul, r[3]_d) (see header comments for locTable)
396  locTable.initialize( 1, 0, 1, 3, num_points );
398 
399  // pass storage directly into locate_points, since we know those arrays are contiguous
400  ErrorCode rval = locate_points( pos, num_points, (EntityHandle*)locTable.vul_wr, locTable.vr_wr, NULL, rel_iter_tol,
401  abs_iter_tol, inside_tol );
402  std::fill( locTable.vi_wr, locTable.vi_wr + num_points, 0 );
403  locTable.set_n( num_points );
404  if( MB_SUCCESS != rval ) return rval;
405 
406  // only call this if I'm the top-level function, since it resets the last time called
408 
409  return MB_SUCCESS;
410 }
411 
413  EntityHandle* ents,
414  double* params,
415  int* is_inside,
416  const double rel_iter_tol,
417  const double abs_iter_tol,
418  const double inside_tol )
419 {
420  bool i_initialized = false;
421  if( !timerInitialized )
422  {
424  timerInitialized = true;
425  i_initialized = true;
426  }
427 
428  assert( !verts.empty() && mbImpl->type_from_handle( *verts.rbegin() ) == MBVERTEX );
429  std::vector< double > pos( 3 * verts.size() );
430  ErrorCode rval = mbImpl->get_coords( verts, &pos[0] );
431  if( MB_SUCCESS != rval ) return rval;
432  rval = locate_points( &pos[0], verts.size(), ents, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
433 
434  // only call this if I'm the top-level function, since it resets the last time called
436 
437  return rval;
438 }
439 
441  int num_points,
442  EntityHandle* ents,
443  double* params,
444  int* is_inside,
445  const double /* rel_iter_tol */,
446  const double abs_iter_tol,
447  const double inside_tol )
448 {
449  bool i_initialized = false;
450  if( !timerInitialized )
451  {
453  timerInitialized = true;
454  i_initialized = true;
455  }
456 
457  /*
458  double tmp_abs_iter_tol = abs_iter_tol;
459  if (rel_iter_tol && !tmp_abs_iter_tol) {
460  // relative epsilon given, translate to absolute epsilon using box dimensions
461  tmp_abs_iter_tol = rel_iter_tol * localBox.diagonal_length();
462  }
463  */
464 
466 
467  ErrorCode rval = MB_SUCCESS;
468  for( int i = 0; i < num_points; i++ )
469  {
470  int i3 = 3 * i;
471  ErrorCode tmp_rval =
472  myTree->point_search( pos + i3, ents[i], abs_iter_tol, inside_tol, NULL, NULL, (CartVect*)( params + i3 ) );
473  if( MB_SUCCESS != tmp_rval )
474  {
475  rval = tmp_rval;
476  continue;
477  }
478 
479  if( debug && !ents[i] )
480  {
481  std::cout << "Point " << i << " not found; point: (" << pos[i3] << "," << pos[i3 + 1] << "," << pos[i3 + 2]
482  << ")" << std::endl;
483  }
484 
485  if( is_inside ) is_inside[i] = ( ents[i] ? true : false );
486  }
487 
488  // only call this if I'm the top-level function, since it resets the last time called
490 
491  return rval;
492 }
493 
494 /* Count the number of located points in locTable
495  * Return the number of entries in locTable that have non-zero entity handles, which
496  * represents the number of points in targetEnts that were inside one element in sourceEnts
497  *
498  */
500 {
501  int num_located = locTable.get_n() - std::count( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
502  if( num_located != (int)locTable.get_n() )
503  {
504  Ulong* nl = std::find( locTable.vul_rd, locTable.vul_rd + locTable.get_n(), 0 );
505  if( nl )
506  {
507  int idx = nl - locTable.vul_rd;
508  if( idx )
509  {
510  }
511  }
512  }
513  return num_located;
514 }
515 
516 /* Count the number of located points in parLocTable
517  * Return the number of entries in parLocTable that have a non-negative index in on a remote
518  * proc in parLocTable, which gives the number of points located in at least one element in a
519  * remote proc's sourceEnts.
520  */
522 {
523  int located = 0;
524  for( unsigned int i = 0; i < parLocTable.get_n(); i++ )
525  if( parLocTable.vi_rd[2 * i] != -1 ) located++;
526  return located;
527 }
528 } // namespace moab