 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "moab/ParallelMergeMesh.hpp"
2 #include "moab/Core.hpp"
3 #include "moab/CartVect.hpp"
4 #include "moab/BoundBox.hpp"
5 #include "moab/Skinner.hpp"
6 #include "moab/MergeMesh.hpp"
7 #include "moab/CN.hpp"
8 #include <cfloat>
9 #include <algorithm>
10
11 #ifdef MOAB_HAVE_MPI
12 #include "moab_mpi.h"
13 #endif
14
15 namespace moab
16 {
17
18 // Constructor
19 /*Get Merge Data and tolerance*/
20 ParallelMergeMesh::ParallelMergeMesh( ParallelComm* pc, const double epsilon ) : myPcomm( pc ), myEps( epsilon )
21 {
22 myMB = pc->get_moab();
23 mySkinEnts.resize( 4 );
24 }
25
26 // Have a wrapper function on the actual merge to avoid memory leaks
27 // Merges elements within a proximity of epsilon
28 ErrorCode ParallelMergeMesh::merge( EntityHandle levelset, bool skip_local_merge, int dim )
29 {
30 ErrorCode rval = PerformMerge( levelset, skip_local_merge, dim );MB_CHK_ERR( rval );
31 CleanUp();
32 return rval;
33 }
34
35 // Perform the merge
36 ErrorCode ParallelMergeMesh::PerformMerge( EntityHandle levelset, bool skip_local_merge, int dim )
37 {
38 // Get the mesh dimension
39 ErrorCode rval;
40 if( dim < 0 )
41 {
42 rval = myMB->get_dimension( dim );MB_CHK_ERR( rval );
43 }
44
45 // Get the local skin elements
46 rval = PopulateMySkinEnts( levelset, dim, skip_local_merge );
47 // If there is only 1 proc, we can return now
48 if( rval != MB_SUCCESS || myPcomm->size() == 1 )
49 {
50 return rval;
51 }
52
53 // Determine the global bounding box
54 double gbox[6];
55 rval = GetGlobalBox( gbox );MB_CHK_ERR( rval );
56
57 /* Assemble The Destination Tuples */
58 // Get a list of tuples which contain (toProc, handle, x,y,z)
59 myTup.initialize( 1, 0, 1, 3, mySkinEnts[0].size() );
60 rval = PopulateMyTup( gbox );MB_CHK_ERR( rval );
61
62 /* Gather-Scatter Tuple
63 -tup comes out as (remoteProc,handle,x,y,z) */
64 myCD.initialize( myPcomm->comm() );
65
66 // 1 represents dynamic tuple, 0 represents index of the processor to send to
67 myCD.gs_transfer( 1, myTup, 0 );
68
69 /* Sort By X,Y,Z
70 -Utilizes a custom quick sort incorporating epsilon*/
71 SortTuplesByReal( myTup, myEps );
72
73 // Initialize another tuple list for matches
74 myMatches.initialize( 2, 0, 2, 0, mySkinEnts[0].size() );
75
76 // ID the matching tuples
77 rval = PopulateMyMatches();MB_CHK_ERR( rval );
78
79 // We can free up the tuple myTup now
80 myTup.reset();
81
82 /*Gather-Scatter Again*/
83 // 1 represents dynamic list, 0 represents proc index to send tuple to
84 myCD.gs_transfer( 1, myMatches, 0 );
85 // We can free up the crystal router now
86 myCD.reset();
87
88 // Sort the matches tuple list
89 SortMyMatches();
90
91 // Tag the shared elements
92 rval = TagSharedElements( dim );MB_CHK_ERR( rval );
93
94 // Free up the matches tuples
95 myMatches.reset();
96 return rval;
97 }
98
99 // Sets mySkinEnts with all of the skin entities on the processor
100 ErrorCode ParallelMergeMesh::PopulateMySkinEnts( const EntityHandle meshset, int dim, bool skip_local_merge )
101 {
102 /*Merge Mesh Locally*/
103 // Get all dim dimensional entities
104 Range ents;
105 ErrorCode rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval );
106
107 if( ents.empty() && dim == 3 )
108 {
109 dim--;
110 rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval ); // maybe dimension 2
111 }
112
113 // Merge Mesh Locally
114 if( !skip_local_merge )
115 {
116 MergeMesh merger( myMB, false );
117 merger.merge_entities( ents, myEps );
118 // We can return if there is only 1 proc
119 if( rval != MB_SUCCESS || myPcomm->size() == 1 )
120 {
121 return rval;
122 }
123
124 // Rebuild the ents range
125 ents.clear();
126 rval = myMB->get_entities_by_dimension( meshset, dim, ents );MB_CHK_ERR( rval );
127 }
128
129 /*Get Skin
130 -Get Range of all dimensional entities
131 -skinEnts[i] is the skin entities of dimension i*/
132 Skinner skinner( myMB );
133 for( int skin_dim = dim; skin_dim >= 0; skin_dim-- )
134 {
135 rval = skinner.find_skin( meshset, ents, skin_dim, mySkinEnts[skin_dim] );MB_CHK_ERR( rval );
136 }
137 return MB_SUCCESS;
138 }
139
140 // Determine the global assembly box
141 ErrorCode ParallelMergeMesh::GetGlobalBox( double* gbox )
142 {
143 ErrorCode rval;
144
145 /*Get Bounding Box*/
146 BoundBox box;
147 if( mySkinEnts[0].size() != 0 )
148 {
149 rval = box.update( *myMB, mySkinEnts[0] );MB_CHK_ERR( rval );
150 }
151
152 // Invert the max
153 box.bMax *= -1;
154
155 /*Communicate to all processors*/
156 MPI_Allreduce( (void*)&box, gbox, 6, MPI_DOUBLE, MPI_MIN, myPcomm->comm() );
157
158 /*Assemble Global Bounding Box*/
159 // Flip the max back
160 for( int i = 3; i < 6; i++ )
161 {
162 gbox[i] *= -1;
163 }
164 return MB_SUCCESS;
165 }
166
167 // Assemble the tuples with their processor destination
168 ErrorCode ParallelMergeMesh::PopulateMyTup( double* gbox )
169 {
170 /*Figure out how do partition the global box*/
171 double lengths[3];
172 int parts[3];
173 ErrorCode rval = PartitionGlobalBox( gbox, lengths, parts );MB_CHK_ERR( rval );
174
175 /* Get Skin Coordinates, Vertices */
176 double* x = new double[mySkinEnts[0].size()];
177 double* y = new double[mySkinEnts[0].size()];
178 double* z = new double[mySkinEnts[0].size()];
179 rval = myMB->get_coords( mySkinEnts[0], x, y, z );
180 if( rval != MB_SUCCESS )
181 {
182 // Prevent Memory Leak
183 delete[] x;
184 delete[] y;
185 delete[] z;
186 return rval;
187 }
188
189 // Initialize variable to be used in the loops
190 std::vector< int > toProcs;
191 int xPart, yPart, zPart, xEps, yEps, zEps, baseProc;
192 unsigned long long tup_i = 0, tup_ul = 0, tup_r = 0, count = 0;
193 // These are boolean to determine if the vertex is on close enough to a given border
194 bool xDup, yDup, zDup;
195 bool canWrite = myTup.get_writeEnabled();
196 if( !canWrite ) myTup.enableWriteAccess();
197 // Go through each vertex
198 for( Range::iterator it = mySkinEnts[0].begin(); it != mySkinEnts[0].end(); ++it )
199 {
200 xDup = false;
201 yDup = false;
202 zDup = false;
203 // Figure out which x,y,z partition the element is in.
204 xPart = static_cast< int >( floor( ( x[count] - gbox[0] ) / lengths[0] ) );
205 xPart = ( xPart < parts[0] ? xPart : parts[0] - 1 ); // Make sure it stays within the bounds
206
207 yPart = static_cast< int >( floor( ( y[count] - gbox[1] ) / lengths[1] ) );
208 yPart = ( yPart < parts[1] ? yPart : parts[1] - 1 ); // Make sure it stays within the bounds
209
210 zPart = static_cast< int >( floor( ( z[count] - gbox[2] ) / lengths[2] ) );
211 zPart = ( zPart < parts[2] ? zPart : parts[2] - 1 ); // Make sure it stays within the bounds
212
213 // Figure out the partition with the addition of Epsilon
214 xEps = static_cast< int >( floor( ( x[count] - gbox[0] + myEps ) / lengths[0] ) );
215 yEps = static_cast< int >( floor( ( y[count] - gbox[1] + myEps ) / lengths[1] ) );
216 zEps = static_cast< int >( floor( ( z[count] - gbox[2] + myEps ) / lengths[2] ) );
217
218 // Figure out if the vertex needs to be sent to multiple procs
219 xDup = ( xPart != xEps && xEps < parts[0] );
220 yDup = ( yPart != yEps && yEps < parts[1] );
221 zDup = ( zPart != zEps && zEps < parts[2] );
222
223 // Add appropriate processors to the vector
224 baseProc = xPart + yPart * parts[0] + zPart * parts[0] * parts[1];
225 toProcs.push_back( baseProc );
226 if( xDup )
227 {
228 toProcs.push_back( baseProc + 1 ); // Get partition to the right
229 }
230 if( yDup )
231 {
232 // Partition up 1
233 toProcs.push_back( baseProc + parts[0] );
234 }
235 if( zDup )
236 {
237 // Partition above 1
238 toProcs.push_back( baseProc + parts[0] * parts[1] );
239 }
240 if( xDup && yDup )
241 {
242 // Partition up 1 and right 1
243 toProcs.push_back( baseProc + parts[0] + 1 );
244 }
245 if( xDup && zDup )
246 {
247 // Partition right 1 and above 1
248 toProcs.push_back( baseProc + parts[0] * parts[1] + 1 );
249 }
250 if( yDup && zDup )
251 {
252 // Partition up 1 and above 1
253 toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] );
254 }
255 if( xDup && yDup && zDup )
256 {
257 // Partition right 1, up 1, and above 1
258 toProcs.push_back( baseProc + parts[0] * parts[1] + parts[0] + 1 );
259 }
260 // Grow the tuple list if necessary
261 while( myTup.get_n() + toProcs.size() >= myTup.get_max() )
262 {
263 myTup.resize( myTup.get_max() ? myTup.get_max() + myTup.get_max() / 2 + 1 : 2 );
264 }
265
266 // Add each proc as a tuple
267 for( std::vector< int >::iterator proc = toProcs.begin(); proc != toProcs.end(); ++proc )
268 {
269 myTup.vi_wr[tup_i++] = *proc;
270 myTup.vul_wr[tup_ul++] = *it;
271 myTup.vr_wr[tup_r++] = x[count];
272 myTup.vr_wr[tup_r++] = y[count];
273 myTup.vr_wr[tup_r++] = z[count];
274 myTup.inc_n();
275 }
276 count++;
277 toProcs.clear();
278 }
279 delete[] x;
280 delete[] y;
281 delete[] z;
282 if( !canWrite ) myTup.disableWriteAccess();
283 return MB_SUCCESS;
284 }
285
286 // Partition the global box by the number of procs
287 ErrorCode ParallelMergeMesh::PartitionGlobalBox( double* gbox, double* lengths, int* parts )
288 {
289 // Determine the length of each side
290 double xLen = gbox[3] - gbox[0];
291 double yLen = gbox[4] - gbox[1];
292 double zLen = gbox[5] - gbox[2];
293 // avoid division by zero for deciding the way to partition the global box
294 // make all sides of the box at least myEps
295 // it can be zero in some cases, so not possible to compare the ratios :) for best division
296 if( xLen < myEps ) xLen = myEps;
297 if( yLen < myEps ) yLen = myEps;
298 if( zLen < myEps ) zLen = myEps;
299 unsigned numProcs = myPcomm->size();
300
301 // Partition sides from the longest to shortest lengths
302 // If x is the longest side
303 if( xLen >= yLen && xLen >= zLen )
304 {
305 parts[0] = PartitionSide( xLen, yLen * zLen, numProcs, true );
306 numProcs /= parts[0];
307 // If y is second longest
308 if( yLen >= zLen )
309 {
310 parts[1] = PartitionSide( yLen, zLen, numProcs, false );
311 parts[2] = numProcs / parts[1];
312 }
313 // If z is the longer
314 else
315 {
316 parts[2] = PartitionSide( zLen, yLen, numProcs, false );
317 parts[1] = numProcs / parts[2];
318 }
319 }
320 // If y is the longest side
321 else if( yLen >= zLen )
322 {
323 parts[1] = PartitionSide( yLen, xLen * zLen, numProcs, true );
324 numProcs /= parts[1];
325 // If x is the second longest
326 if( xLen >= zLen )
327 {
328 parts[0] = PartitionSide( xLen, zLen, numProcs, false );
329 parts[2] = numProcs / parts[0];
330 }
331 // If z is the second longest
332 else
333 {
334 parts[2] = PartitionSide( zLen, xLen, numProcs, false );
335 parts[0] = numProcs / parts[2];
336 }
337 }
338 // If z is the longest side
339 else
340 {
341 parts[2] = PartitionSide( zLen, xLen * yLen, numProcs, true );
342 numProcs /= parts[2];
343 // If x is the second longest
344 if( xLen >= yLen )
345 {
346 parts[0] = PartitionSide( xLen, yLen, numProcs, false );
347 parts[1] = numProcs / parts[0];
348 }
349 // If y is the second longest
350 else
351 {
352 parts[1] = PartitionSide( yLen, xLen, numProcs, false );
353 parts[0] = numProcs / parts[1];
354 }
355 }
356
357 // Divide up each side to give the lengths
358 lengths[0] = xLen / (double)parts[0];
359 lengths[1] = yLen / (double)parts[1];
360 lengths[2] = zLen / (double)parts[2];
361 return MB_SUCCESS;
362 }
363
364 // Partition a side based on the length ratios
365 int ParallelMergeMesh::PartitionSide( double sideLen, double restLen, unsigned numProcs, bool altRatio )
366 {
367 // If theres only 1 processor, then just return 1
368 if( numProcs == 1 )
369 {
370 return 1;
371 }
372 // Initialize with the ratio of 1 proc
373 double ratio = -DBL_MAX;
374 unsigned factor = 1;
375 // We need to be able to save the last ratio and factor (for comparison)
376 double oldRatio = ratio;
377 double oldFactor = 1;
378
379 // This is the ratio were shooting for
380 double goalRatio = sideLen / restLen;
381
382 // Calculate the divisor and numerator power
383 // This avoid if statements in the loop and is useful since both calculations are similar
384 double divisor, p;
385 if( altRatio )
386 {
387 divisor = (double)numProcs * sideLen;
388 p = 3;
389 }
390 else
391 {
392 divisor = (double)numProcs;
393 p = 2;
394 }
395
396 // Find each possible factor
397 for( unsigned i = 2; i <= numProcs / 2; i++ )
398 {
399 // If it is a factor...
400 if( numProcs % i == 0 )
401 {
402 // We need to save the past factor
403 oldRatio = ratio;
404 oldFactor = factor;
405 // There are 2 different ways to calculate the ratio:
406 // Comparing 1 side to 2 sides: (i*i*i)/(numProcs*x)
407 // Justification: We have a ratio x:y:z (side Lengths) == a:b:c (procs). So a=kx,
408 // b=ky, c=kz. Also, abc=n (numProcs) => bc = n/a. Also, a=kx => k=a/x => 1/k=x/a And so
409 // x/(yz) == (kx)/(kyz) == (kx)/(kykz(1/k)) == a/(bc(x/a)) == a/((n/a)(x/a)) == a^3/(nx).
410 // Comparing 1 side to 1 side: (i*i)/numprocs
411 // Justification: i/(n/i) == i^2/n
412 ratio = pow( (double)i, p ) / divisor;
413 factor = i;
414 // Once we have passed the goal ratio, we can break since we'll only move away from the
415 // goal ratio
416 if( ratio >= goalRatio )
417 {
418 break;
419 }
420 }
421 }
422 // If we haven't reached the goal ratio yet, check out factor = numProcs
423 if( ratio < goalRatio )
424 {
425 oldRatio = ratio;
426 oldFactor = factor;
427 factor = numProcs;
428 ratio = pow( (double)numProcs, p ) / divisor;
429 }
430
431 // Figure out if our oldRatio is better than ratio
432 if( fabs( ratio - goalRatio ) > fabs( oldRatio - goalRatio ) )
433 {
434 factor = oldFactor;
435 }
436 // Return our findings
437 return factor;
438 }
439
440 // Id the tuples that are matching
441 ErrorCode ParallelMergeMesh::PopulateMyMatches()
442 {
443 // Counters for accessing tuples more efficiently
444 unsigned long i = 0, mat_i = 0, mat_ul = 0, j = 0, tup_r = 0;
445 double eps2 = myEps * myEps;
446
447 uint tup_mi, tup_ml, tup_mul, tup_mr;
448 myTup.getTupleSize( tup_mi, tup_ml, tup_mul, tup_mr );
449
450 bool canWrite = myMatches.get_writeEnabled();
451 if( !canWrite ) myMatches.enableWriteAccess();
452
453 while( ( i + 1 ) < myTup.get_n() )
454 {
455 // Proximity Comparison
456 double xi = myTup.vr_rd[tup_r], yi = myTup.vr_rd[tup_r + 1], zi = myTup.vr_rd[tup_r + 2];
457
458 bool done = false;
459 while( !done )
460 {
461 j++;
462 tup_r += tup_mr;
463 if( j >= myTup.get_n() )
464 {
465 break;
466 }
467 CartVect cv( myTup.vr_rd[tup_r] - xi, myTup.vr_rd[tup_r + 1] - yi, myTup.vr_rd[tup_r + 2] - zi );
468 if( cv.length_squared() > eps2 )
469 {
470 done = true;
471 }
472 }
473 // Allocate the tuple list before adding matches
474 while( myMatches.get_n() + ( j - i ) * ( j - i - 1 ) >= myMatches.get_max() )
475 {
476 myMatches.resize( myMatches.get_max() ? myMatches.get_max() + myMatches.get_max() / 2 + 1 : 2 );
477 }
478
479 // We now know that tuples [i to j) exclusive match.
480 // If n tuples match, n*(n-1) match tuples will be made
481 // tuples are of the form (proc1,proc2,handle1,handle2)
482 if( i + 1 < j )
483 {
484 int kproc = i * tup_mi;
485 unsigned long khand = i * tup_mul;
486 for( unsigned long k = i; k < j; k++ )
487 {
488 int lproc = kproc + tup_mi;
489 unsigned long lhand = khand + tup_mul;
490 for( unsigned long l = k + 1; l < j; l++ )
491 {
492 myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc1
493 myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc2
494 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle1
495 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle2
496 myMatches.inc_n();
497
498 myMatches.vi_wr[mat_i++] = myTup.vi_rd[lproc]; // proc1
499 myMatches.vi_wr[mat_i++] = myTup.vi_rd[kproc]; // proc2
500 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[lhand]; // handle1
501 myMatches.vul_wr[mat_ul++] = myTup.vul_rd[khand]; // handle2
502 myMatches.inc_n();
503 lproc += tup_mi;
504 lhand += tup_mul;
505 }
506 kproc += tup_mi;
507 khand += tup_mul;
508 } // End for(int k...
509 }
510 i = j;
511 } // End while(i+1<tup.n)
512
513 if( !canWrite ) myMatches.disableWriteAccess();
514 return MB_SUCCESS;
515 }
516
517 // Sort the matching tuples so that vertices can be tagged accurately
518 ErrorCode ParallelMergeMesh::SortMyMatches()
519 {
520 TupleList::buffer buf( mySkinEnts[0].size() );
521 // Sorts are necessary to check for doubles
522 // Sort by remote handle
523 myMatches.sort( 3, &buf );
524 // Sort by matching proc
525 myMatches.sort( 1, &buf );
526 // Sort by local handle
527 myMatches.sort( 2, &buf );
528 buf.reset();
529 return MB_SUCCESS;
530 }
531
532 // Tag the shared elements using existing PComm functionality
533 ErrorCode ParallelMergeMesh::TagSharedElements( int dim )
534 {
535 // Manipulate the matches list to tag vertices and entities
536 // Set up proc ents
537 Range proc_ents;
538 ErrorCode rval;
539
540 // get the entities in the partition sets
541 for( Range::iterator rit = myPcomm->partitionSets.begin(); rit != myPcomm->partitionSets.end(); ++rit )
542 {
543 Range tmp_ents;
544 rval = myMB->get_entities_by_handle( *rit, tmp_ents, true );
545 if( MB_SUCCESS != rval )
546 {
547 return rval;
548 }
549 proc_ents.merge( tmp_ents );
550 }
551 if( myMB->dimension_from_handle( *proc_ents.rbegin() ) != myMB->dimension_from_handle( *proc_ents.begin() ) )
552 {
553 Range::iterator lower = proc_ents.lower_bound( CN::TypeDimensionMap[0].first ),
554 upper = proc_ents.upper_bound( CN::TypeDimensionMap[dim - 1].second );
555 proc_ents.erase( lower, upper );
556 }
557
558 // This vector doesn't appear to be used but its in resolve_shared_ents
559 int maxp = -1;
560 std::vector< int > sharing_procs( MAX_SHARING_PROCS );
561 std::fill( sharing_procs.begin(), sharing_procs.end(), maxp );
562
563 // get ents shared by 1 or n procs
564 std::map< std::vector< int >, std::vector< EntityHandle > > proc_nranges;
565 Range proc_verts;
566 rval = myMB->get_adjacencies( proc_ents, 0, false, proc_verts, Interface::UNION );
567 if( rval != MB_SUCCESS )
568 {
569 return rval;
570 }
571
572 rval = myPcomm->tag_shared_verts( myMatches, proc_nranges, proc_verts );
573 if( rval != MB_SUCCESS )
574 {
575 return rval;
576 }
577
578 // get entities shared by 1 or n procs
579 rval = myPcomm->get_proc_nvecs( dim, dim - 1, &mySkinEnts[0], proc_nranges );
580 if( rval != MB_SUCCESS )
581 {
582 return rval;
583 }
584
585 // create the sets for each interface; store them as tags on
586 // the interface instance
587 Range iface_sets;
588 rval = myPcomm->create_interface_sets( proc_nranges );
589 if( rval != MB_SUCCESS )
590 {
591 return rval;
592 }
593 // establish comm procs and buffers for them
594 std::set< unsigned int > procs;
595 rval = myPcomm->get_interface_procs( procs, true );
596 if( rval != MB_SUCCESS )
597 {
598 return rval;
599 }
600
601 // resolve shared entity remote handles; implemented in ghost cell exchange
602 // code because it's so similar
603 rval = myPcomm->exchange_ghost_cells( -1, -1, 0, true, true );
604 if( rval != MB_SUCCESS )
605 {
606 return rval;
607 }
608 // now build parent/child links for interface sets
609 rval = myPcomm->create_iface_pc_links();
610 return rval;
611 }
612
613 // Make sure to free up any allocated data
614 // Need to avoid a double free
615 void ParallelMergeMesh::CleanUp()
616 {
617 // The reset operation is now safe and avoids a double free()
618 myMatches.reset();
619 myTup.reset();
620 myCD.reset();
621 }
622
623 // Simple quick sort to real
624 void ParallelMergeMesh::SortTuplesByReal( TupleList& tup, double eps )
625 {
626 bool canWrite = tup.get_writeEnabled();
627 if( !canWrite ) tup.enableWriteAccess();
628
629 uint mi, ml, mul, mr;
630 tup.getTupleSize( mi, ml, mul, mr );
631 PerformRealSort( tup, 0, tup.get_n(), eps, mr );
632
633 if( !canWrite ) tup.disableWriteAccess();
634 }
635
636 // Swap around tuples
637 void ParallelMergeMesh::SwapTuples( TupleList& tup, unsigned long a, unsigned long b )
638 {
639 if( a == b ) return;
640
641 uint mi, ml, mul, mr;
642 tup.getTupleSize( mi, ml, mul, mr );
643
644 // Swap mi
645 unsigned long a_val = a * mi, b_val = b * mi;
646 for( unsigned long i = 0; i < mi; i++ )
647 {
648 sint t = tup.vi_rd[a_val];
649 tup.vi_wr[a_val] = tup.vi_rd[b_val];
650 tup.vi_wr[b_val] = t;
651 a_val++;
652 b_val++;
653 }
654 // Swap ml
655 a_val = a * ml;
656 b_val = b * ml;
657 for( unsigned long i = 0; i < ml; i++ )
658 {
659 slong t = tup.vl_rd[a_val];
660 tup.vl_wr[a_val] = tup.vl_rd[b_val];
661 tup.vl_wr[b_val] = t;
662 a_val++;
663 b_val++;
664 }
665 // Swap mul
666 a_val = a * mul;
667 b_val = b * mul;
668 for( unsigned long i = 0; i < mul; i++ )
669 {
670 Ulong t = tup.vul_rd[a_val];
671 tup.vul_wr[a_val] = tup.vul_rd[b_val];
672 tup.vul_wr[b_val] = t;
673 a_val++;
674 b_val++;
675 }
676 // Swap mr
677 a_val = a * mr;
678 b_val = b * mr;
679 for( unsigned long i = 0; i < mr; i++ )
680 {
681 realType t = tup.vr_rd[a_val];
682 tup.vr_wr[a_val] = tup.vr_rd[b_val];
683 tup.vr_wr[b_val] = t;
684 a_val++;
685 b_val++;
686 }
687 }
688
689 // Perform the sorting of a tuple by real
690 // To sort an entire tuple_list, call (tup,0,tup.n,epsilon)
691 void ParallelMergeMesh::PerformRealSort( TupleList& tup,
692 unsigned long left,
693 unsigned long right,
694 double eps,
695 uint tup_mr )
696 {
697 // If list size is only 1 or 0 return
698 if( left + 1 >= right )
699 {
700 return;
701 }
702 unsigned long swap = left, tup_l = left * tup_mr, tup_t = tup_l + tup_mr;
703
704 // Swap the median with the left position for a (hopefully) better split
705 SwapTuples( tup, left, ( left + right ) / 2 );
706
707 // Partition the data
708 for( unsigned long t = left + 1; t < right; t++ )
709 {
710 // If the left value(pivot) is greater than t_val, swap it into swap
711 if( TupleGreaterThan( tup, tup_l, tup_t, eps, tup_mr ) )
712 {
713 swap++;
714 SwapTuples( tup, swap, t );
715 }
716 tup_t += tup_mr;
717 }
718
719 // Swap so that position swap is in the correct position
720 SwapTuples( tup, left, swap );
721
722 // Sort left and right of swap
723 PerformRealSort( tup, left, swap, eps, tup_mr );
724 PerformRealSort( tup, swap + 1, right, eps, tup_mr );
725 }
726
727 // Note, this takes the actual tup.vr[] index (aka i*tup.mr)
728 bool ParallelMergeMesh::TupleGreaterThan( TupleList& tup,
729 unsigned long vrI,
730 unsigned long vrJ,
731 double eps,
732 uint tup_mr )
733 {
734 unsigned check = 0;
735 while( check < tup_mr )
736 {
737 // If the values are the same
738 if( fabs( tup.vr_rd[vrI + check] - tup.vr_rd[vrJ + check] ) <= eps )
739 {
740 check++;
741 continue;
742 }
743 // If I greater than J
744 else if( tup.vr_rd[vrI + check] > tup.vr_rd[vrJ + check] )
745 {
746 return true;
747 }
748 // If J greater than I
749 else
750 {
751 return false;
752 }
753 }
754 // All Values are the same
755 return false;
756 }
757
758 } // End namespace moab
1.9.1.