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